1

我对 R 中的以下函数有疑问：

test <- function(alpha, beta, n){
  result <- exp(lgamma(alpha) + lgamma(n + beta) - lgamma(alpha + beta + n) - (lgamma(alpha) + lgamma(beta) - lgamma(alpha + beta)))
  return(result)
}

现在，如果您插入以下值：

betabinom(-0.03292708, -0.3336882, 10)

它应该失败并导致NaN. 那是因为如果我们在 Excel 中实现精确的函数，我们会得到一个不是数字的结果。Excel 中的实现很简单，因为 J32 是alpha、K32beta和 L32的单元格N。生成的单元格的实现如下所示：

=EXP(GAMMALN(J32)+GAMMALN(L32+K32)-GAMMALN(J32+K32+L32)-(GAMMALN(J32)+GAMMALN(K32)-GAMMALN(J32+K32)))

所以这似乎给出了正确的答案，因为该函数仅针对大于零的 alpha 和 beta 以及大于或等于零的 n 定义。因此我想知道这里发生了什么？我还尝试了 Rmpf 包来提高数值精度，但这似乎没有任何作用。

谢谢

4

1 回答 1

2

tl;dr log(gamma(x)) 的定义比您想象的或 Excel 认为的更普遍。alpha如果您希望您的函数不接受and的负值beta或返回NaN，只需手动测试并返回适当的值 ( if (alpha<0 || beta<0) return(NaN))。

这不是数值精度问题，而是定义问题。Gamma 函数是为负实数值定义的：?lgamma说：

伽马函数由（Abramowitz 和 Stegun 部分 6.1.1，第 255 页）定义

Gamma(x) = 积分_0^Inf t^(x-1) exp(-t) dt

对于除零和负整数以外的所有实数“x”（当返回“NaN”时）。

此外，参考lgamma...

...以及伽马函数绝对值的自然对数...

（重点是原文）

curve(lgamma(x),-1,1)

在此处输入图像描述

gamma(-0.1)          ## -10.68629
log(gamma(-0.1)+0i)  ## 2.368961+3.141593i
log(abs(gamma(-0.1)) ## 2.368961
lgamma(-0.1)         ## 2.368961

Wolfram Alpha 同意第二个计算。

于 2015-08-06T12:49:55.000 回答