r - R中的多精度伽玛函数

Question

我需要计算 R 中涉及每个元素的 gamma 函数的总和。当 gamma 的参数增加时，我得到 NaN，我怀疑问题是 gamma 评估的数值问题。我已经阅读了 Rmpfr 和 gmp 的文档，但我只发现整数的阶乘。我也在这里发布代码，也许您对错误的来源有更好的了解。

 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 是一个总和为 1 的正实数数组，a1 和 b1 是正值向量。当 a1（和 a1mat）变长并具有高值（~150）时，A 变为 NaN。

Answer 1

尝试在日志空间中工作：

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