使用递归实现 log-gamma

Question

我正在尝试使用递归计算 lgamma。但它没有按预期工作并给出 Nan 和 -Inf。这可能是什么原因？

sum = 0
log_gamma_recursive <- function(n) {
  if(n == 1){
    return(1)
  }else{
    sum = sum + log(log_gamma_recursive(n-1))
    print(sum)
    return(sum)
  }
}#log 4  + log 3 + log 2 + log 1
log_gamma_recursive(5)

Answer 1

代码有三处错误：

sum 不应该是一个全局变量，因为连续的 "external" 调用将 return 一个累加值
您将 return 值的 log 添加到 sum。这将给出一系列 nested 日志：
- 你得到log 4 + log(log 3 + log(log 2 + log (log 1))))，这显然是错误的
- log(log 1) = log 0 从数学上讲它趋向于负无穷大，但是库可能会也可能不会将其归类为有效输入
log_gamma(1) 是 0，不是 1

考虑到这一点，试试这个：

log_gamma_recursive <- function(n) {
  if (n <= 1) {
    return (0)
  }else{
    sum = (log(n-1) + log_gamma_recursive(n-1))
    print(sum)
    return(sum)
  }
}

Answer 2

由于R可以进行向量化计算，所以对结果向量进行log和求和是安全的。

 ll=function(x)sum(log(1:(x-1)))

这样做有什么好处？这可以计算 log_gamma_recursive 无法计算的值。 Recall log_gamma_recursive 是一个嵌套函数，因此在某些时候会失败，具体取决于机器容量。例如：

log_gamma_recursive(3476)
[1] 24862.89
ll(3476)
[1] 24862.89

在这台机器上，log_gamma_recursive 对任何大于 3476 的数字都不起作用：

log_gamma_recursive(3477)
Error: C stack usage  7970456 is too close to the limit

而另一方面，ll 函数有效：

ll(3477)
[1] 24871.04

同时，ll函数远比log_recursive函数快很多：

system.time(for(i in 1:10000)ll(3000))
   user  system elapsed 
  0.313   0.036   0.350 
system.time(for(i in 1:10000)log_gamma_recursive(3000))
   user  system elapsed 
 20.569   0.100  20.794

这些单位是秒！！

使用递归实现 log-gamma

Implenting log-gamma using recursion

r

recursion

gamma