Gamma函数负输入致R中分布均值计算失败，求解决方案

我在计算某一分布的均值，Python实现可顺利运行且结果与参考文献一致，但R版本无法得到正确结果。经排查，问题源于Gamma函数接收了负输入——b-i在多数情况下为负值。手头该分布相关资源多为R语言版本，无需Python实现，仅需解决R中的问题，请问是否有可用的R包或修改方法？

Python代码实现 #

import math
import numpy as np
from scipy.special import gamma
def kum_mom(lambd, b, beta, r):
    limit = 100
    i = 0
    sum = 0
    while i <= limit:
        temp=(((-1)**i)*lambd*b*beta*gamma(b)*gamma(1-(r/beta)))/(math.factorial(i)*beta*((lambd*(i+1))**(1-(r/beta)))*gamma(b-i))
        sum = sum + temp
        i = i + 1
    return sum

kum_mom(1, 1, 3, 1)
# 输出：1.3541179394264005

R代码实现及报错情况 #

library(pracma)
kum_mom <- function(lambda, b, beta, r) {
  limit <- 100
  i <- 0
  sum<- 0
  while (i <= limit) {
    temp <- ((-1)^i * lambda * b * beta * gamma(b) * gamma(1 - (r / beta))) / (factorial(i) * beta * ((lambda * (i + 1))^(1 - (r / beta))) * gamma(b - i))
    sum <- sum + temp
    i <- i + 1
  }
  return(sum)
}

kum_mom(1, 1, 3, 1)
# 输出：[1] NaN
# 警告信息：
# 1: In gamma(b - i) : NaNs produced
# 2: In gamma(b - i) : NaNs produced
# ...（共50+条警告）

解决方案：利用Gamma函数欧拉反射公式处理负输入 #

R原生的gamma()函数（包括pracma包的版本）在输入为负数时会返回NaN，但我们可以用欧拉反射公式将负输入的Gamma函数转换为正输入的计算：
$$\Gamma(z)\Gamma(1-z) = \frac{\pi}{\sin(\pi z)}$$
变形后得到：
$$\Gamma(z) = \frac{\pi}{\sin(\pi z)\cdot\Gamma(1-z)}$$

基于此，我们可以在R中自定义一个支持负输入的Gamma函数，替换原代码中的gamma()调用：

修改后的R代码 #

library(pracma)

# 自定义支持负输入的Gamma函数
gamma_extended <- function(z) {
  if (z > 0) {
    return(gamma(z))
  } else {
    # 应用欧拉反射公式
    return(pi / (sin(pi * z) * gamma(1 - z)))
  }
}

kum_mom <- function(lambda, b, beta, r) {
  limit <- 100
  i <- 0
  sum <- 0
  while (i <= limit) {
    temp <- ((-1)^i * lambda * b * beta * gamma_extended(b) * gamma_extended(1 - (r / beta))) / 
      (factorial(i) * beta * ((lambda * (i + 1))^(1 - (r / beta))) * gamma_extended(b - i))
    sum <- sum + temp
    i <- i + 1
  }
  return(sum)
}

# 测试运行
kum_mom(1, 1, 3, 1)
# 输出：[1] 1.354118

这个修改后的代码会正确处理b-i为负值的情况，计算结果与Python版本一致。

内容的提问来源于stack exchange，提问作者César Augusto