R语言自定义超几何分布函数调用返回NaN结果问题求解

问题1：自定义组合数公式语法错误 #

你写的comb函数存在运算优先级逻辑错误，组合数的标准公式为 $C(n,r) = \frac{n!}{r! \times (n-r)!}$，而你的代码中乘法与除法为同级左结合运算，实际执行逻辑为 $\frac{n!}{r!} \times (n-r)!$，和正确公式完全不符。

问题2：大数阶乘溢出 #

即使修正公式语法，直接计算1100这类大数的阶乘也会超出R双精度浮点数的存储上限（factorial(170)就已经返回Inf），最终出现Inf/Inf的计算结果，返回NaN。 #

修复方案 #

方案1：直接使用R内置优化的组合数函数（推荐） #

R自带的choose()函数底层做了数值稳定性优化，不会出现大数溢出问题，直接替换自定义comb函数即可：

hyper.formula <- function(q,m,n,k) {
  choose(m,q)*choose(n,k-q)/choose(m+n,k)
}

调用测试hyper.formula(7,1100,900,10)，可正常返回概率结果约为0.1523。

方案2：修正自定义comb函数并优化数值稳定性 #

如果需要保留自定义实现，可修改公式语法并通过对数运算避免直接计算大数阶乘：

# 修正后的comb函数
comb <- function(n,r) {
  # 用对数阶乘计算再转指数，避免大数溢出
  exp(lfactorial(n) - (lfactorial(r) + lfactorial(n-r)))
}
# 超几何函数不需要修改
hyper.formula <- function(q,m,n,k) {
  comb(m,q)*comb(n,k-q)/comb(m+n,k)
}

内容的提问来源于stack exchange，提问作者moth