R中的多精度伽玛函数

Question

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

 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

尝试在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

R中的多精度伽玛函数

Multiple precision Gamma function in R

r

numeric

gamma-function