R中的Eval非常大的阶乘

Question

我需要整合一个在其表达式中具有阶乘的函数。但是，如果您尝试计算阶乘，当 n > 170 时，R 将返回 Inf。

我发现很多包可以让你计算非常大的数字，但是，它们总是从我无法集成的类中返回一个对象。积分的最终结果总是一个小数。

这是我的代码：

integrand <- function(n, i, x) {
    (factorial(n) / (factorial(i - 1) * factorial(n - i))) *
        x^(i - 1) * (1 - x)^(n - i)
}

forder <- function(Fx, x, n, i, ...) {
    lower <- sapply(x - 1, Fx, ...)
    upper <- sapply(x, Fx, ...)
     integrate(integrand,
               lower = lower, upper = upper, n = n, i = i,
               stop.on.error = FALSE)$value
}
forder <- Vectorize(forder, "x")

##------------------------------------------------------------------------------
## Some example
y <- sort(rpois(100, 1))

## Works fine
forder(ppois, y, 170, 10, lambda = 1)

## Does not work
forder(ppois, y, 171, 10, lambda = 1)
##------------------------------------------------------------------------------

Answer 1

正如我在评论中所说，您可以替换(factorial(n) / (factorial(i - 1) * factorial(n - i)))为i*choose(n, i). 这两个量是相等的，但choose(n,i)允许更高的 值n。

或者您可以使用该pbeta函数而不是进行数值积分：

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  i*choose(n, i) * (pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)) * beta(i, n-i+1)
}

更好的是，使用对数：

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  lg <- log(i) + lchoose(n, i) + 
    log(pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)) + lbeta(i, n-i+1)
  exp(lg)
}

编辑

我没有注意到这种简化：i*choose(n, i) * beta(i, n-i+1) = 1. 所以你可以简单地做：

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)
}