R中的数值精度问题？

Question

我对 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 的单元格，K32 beta 和 L32 是 N 的单元格。结果单元格的实现如下：

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

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

谢谢

Answer 1

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

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

The gamma function is defined by (Abramowitz and Stegun section 6.1.1, page 255)

Gamma(x) = integral_0^Inf t^(x-1) exp(-t) dt

for all real ‘x’ except zero and negative integers (when ‘NaN’ is returned).

另外参考lgamma ...

... and the natural logarithm of the absolute value of the gamma function ...

（原文着重）

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 agrees with second calculation.

R中的数值精度问题？

Numerical precision problems in R?

r

gamma-function