在R 中使用较大的有限上界时，我发现 integrate() 的行为与 Inf 不一致。

考虑下列积分：

integrate(\(x) x^2 * exp(-x), 0, Inf)
# 2 with absolute error < 7.1e-05

这是正确的（真实值为2）。

然而，当我把 Inf 替换成一个较大的有限数时：

integrate(\(x) x^2 * exp(-x), 0, 10000)
# 2 with absolute error < 0.00011

integrate(\(x) x^2 * exp(-x), 0, 100000)
# 2.429968e-41 with absolute error < 4.8e-41

integrate(\(x) x^2 * exp(-x), 0, 1000000)
# 0 with absolute error < 0

在某些点，结果几乎等于零，这显然不正确。

有趣的是，再次减小上界后又能工作：

integrate(\(x) x^2 * exp(-x), 0, 1000)
# 2 with absolute error < 1.6e-06

问题：

1. 为什么对于非常大的有限上界，integrate() 会失败，但对 Inf 却能正确工作？
2. 这是由于数值下溢、积分算法，还是其他原因？
3. 有没有可靠处理这类积分的推荐方法？

关于 integrate() 的内部行为的一些见解，请参阅 这里。

解决方案

help("integrate") 已经提醒我们不要那样做：

对于无穷区间的积分，应显式进行，而不是仅仅把一个很大的数作为端点。

正如你所展示的，这是出于正当的理由：

library(stats)

integrate(dnorm, 0, .Machine$integer.max)
#> 0 with absolute error < 0

integrate(dnorm, 0, Inf)
#> 0.5 with absolute error < 4.7e-05

嗯，你可以使用任何比 .Machine$integer.max*1e299 更大的数，但这有何必要？！

看看 源代码，我们可以看到发生了什么：

if(is.finite(lower) && is.finite(upper)) {
    wk <- .External(C_call_dqags, ...)  # Finite interval
} else { # indefinite integral
    wk <- .External(C_call_dqagi, ...)  # Infinite interval
}

• 有限界：C_call_dqags（QUADPACK的 dqags）
• 任何无限界：C_call_dqagi（QUADPACK的 dqagi）

我们可以进一步深入查看这些函数及其背后的数学，但我并不是这方面的数学专家。不过文档给了我们更多信息：

如果一个或两个极限是无限的，无限区间被映射到一个有限区间。

对于有限区间，使用全局自适应区间划分，并结合Wynn的 Epsilon外推算法，基本步骤是高斯–克龙德求积。

...

当对无限区间进行积分时，应显式处理，而不是仅仅用一个很大的数作为端点。这会增加得到正确答案的概率——对在无限区间上的积分有限的函数，在该区间的大部分区域应接近于零。

对于在有限点集上的取值要公正地反映该函数在其他地方的行为，函数需要有良好的性质，例如可微，除了可能只有少量跳跃或可积奇点之外。

关于你的第三个问题：

有哪些可靠的方法来处理这类积分？

我不确定是否还有更好的方法，但回想大学时光，我们可以把x 替换为u：

library(stats)

f_original <- \(x) x^2 * exp(-x)

## u = exp(-x)
## x = -log(u)
## dx = -1/u du
## then x = 0 -> u = 1
##      x = Inf -> u = 0
f_transformed <- \(u) {
  x = -log(u)
  x^2 * exp(-x) * (-1/exp(-x))
}

## these give the "same" answer
integrate(f_original, 0, Inf)
#> 2 with absolute error < 7.1e-05
integrate(f_transformed, 1, 0)
#> 2 with absolute error < 0.00017

integrate(f_transformed, exp(-0), exp(-1000000))
#> 2 with absolute error < 0.00017
integrate(f_original, 0, 1000000) ## wrong value
#> 0 with absolute error < 0

创建于2026-03-31，使用 reprex v2.1.1