评估 R 中的非常大的阶乘

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)
}

Answer 2

将integrand更改为使用对数，两个调用都有效。我还发布了函数@StéphaneLaurent使用choose/lchoose的想法和pbeta/beta函数。

integrand <- function(n, i, x) {
y <- lfactorial(n) - lfactorial(i - 1) - lfactorial(n - i) +
(i - 1)*log(x) + log(1 - x)*(n - i)
exp(y)
}
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")
integrandSL <- function(n, i, x) {
y <- log(i) + lchoose(n, i) + (i - 1)*log(x) + log(1 - x)*(n - i)
exp(y)
}
forderSL <- function(Fx, x, n, i, ...) {
lower <- sapply(x - 1, Fx, ...)
upper <- sapply(x, Fx, ...)
integrate(integrandSL,
lower = lower, upper = upper, n = n, i = i,
stop.on.error = FALSE)$value
}
forderSL <- Vectorize(forderSL, "x")
forderSL2 <- 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)
}

现在是测试。所有结果都是all.equal。

##-------------------------------------------------
set.seed(1234)
y <- sort(rpois(100, 1))
res_170 <- forder(ppois, y, 170, 10, lambda = 1)
res_171 <- forder(ppois, y, 171, 10, lambda = 1)
resSL_170 <- forderSL(ppois, y, 170, 10, lambda = 1)
resSL_171 <- forderSL(ppois, y, 171, 10, lambda = 1)
resSL2_170 <- forderSL2(ppois, y, 170, 10, lambda = 1)
resSL2_171 <- forderSL2(ppois, y, 171, 10, lambda = 1)
all.equal(res_170, resSL_170)    # TRUE
all.equal(res_170, resSL2_170)   # TRUE
all.equal(res_171, resSL_171)    # TRUE
all.equal(res_171, resSL2_171)   # TRUE