0

I need to compute a sum in R that involves the gamma function for each element. When the arguments of the gamma increase I get NaN as a result, and I suspect that the problem is numerical with the evaluation of the gamma. I already read the documentation of Rmpfr and gmp but I only found factiorial for integers. I also post the code here, maybe you have a better idea about the source of the error.

 V1<-w1*(b1^a1)/gamma(a1)
 VV1<-outer(V1,V1)
 a1mat<-outer(a1, a1-1, FUN="+")
 b1mat<-outer(b1, b1, FUN="+")
 A<-sum(VV1*gamma(a1mat)/(b1mat^a1mat))

w1 is an array of positive real numbers that sums up to 1, a1 and b1 are vectors of positive values. A becames NaN when a1 (and a1mat) becomes long and with high values (~150).

Phil
  • 7,287
  • 3
  • 36
  • 66
MCMP
  • 3
  • 2
  • An array of positive integers that sums to 1? Can you give an example? – Bill O'Brien Feb 14 '22 at 15:35
  • example: w1=(0.2, 0.3, 0.2, 0.1, 0.1, 0.1), a1=(3, 4, 5, 6, 7, 8) b1=(4,4,4,4,4,4). In this situation works. It fails when w1 becames long (and so small numers close to 0), a1 becames longer (around 150). Oh and of course len(w1)=len(a1)=len(b1) @BillO'Brien – MCMP Feb 14 '22 at 15:42
  • 1
    I'm not sure what the context is here, but if you're needing to approximate an integral, maybe rephrasing this in terms of incomplete gamma functions, or the cdf of a relevant distribution, can avoid the explicit calculation of gamma(x) for large x. Also, maybe you can compute log(gamma(x)) instead; log(gamma(x)) is often a separate function in math libraries. Also, maybe try to compute the limit of the summand for large arguments and approximate tail of summation = limit times number of omitted terms. I use Maxima (https://maxima.sourceforge.io) for stuff like that, maybe try Sympy too. – Robert Dodier Feb 14 '22 at 16:38
  • Can you provide an example that you are trying to calculate where it fails? – jblood94 Feb 14 '22 at 17:45

1 Answers1

0

Try working in log-space:

w1 <- c(0.2, 0.3, 0.2, 0.1, 0.1, 0.1)
a1 <- 3:8
b1 <- rep(4, 6)
a1mat <- outer(a1, a1 - 1, "+")
b1mat <- outer(b1, b1, "+")

# working in log-space
logV1 <- log(w1) + a1*log(b1) - lgamma(a1)
logVV1 <- outer(logV1, logV1, "+")
sum(exp(logVV1 + lgamma(a1mat) - a1mat*log(b1mat)))
#> [1] 0.4614941

# compared to original
V1 <- w1*(b1^a1)/gamma(a1)
VV1 <- outer(V1, V1)
sum(VV1*gamma(a1mat)/(b1mat^a1mat))
#> [1] 0.4614941
jblood94
  • 10,340
  • 1
  • 10
  • 15