R中的数值精度问题？答案

【问题标题】：Numerical precision problems in R?R中的数值精度问题？
【发布时间】：2015-08-06 10:59:39
【问题描述】：

我对 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 包来提高数值精度，但这似乎没有任何作用。

谢谢

【问题讨论】：

这个问题没有任何意义。你可能想改写它。

标签： r gamma-function

【解决方案1】：

tl;dr log(gamma(x)) 的定义比您想象的或 Excel 认为的更普遍。如果您希望您的函数不接受 alpha 和 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 agrees with second calculation.

【讨论】：