python - python 中的函数 "integrate.quad"和 R 中的 "integral",'integrate' 给出错误的结果

Question

我需要评估几个积分，我正在使用正常 (0,1) 密度来测试。

在 python 中

import scipy.integrate as integrate
import scipy.stats
import numpy as np

def integrand(x):
    return scipy.stats.norm(0, 1).pdf(x)

print integrate.quad(integrand, -float('inf'), 0)
print integrate.quad(integrand,-np.inf,100)

(0.4999999999999999, 5.089095674629994e-09)

(0.0, 0.0)

令我非常困惑的是，计算机正确计算了范围 (-inf,0) 的积分，但完全错过了 (-inf,100)(应该接近 1)。因此，我在 R 中尝试了以下内容

integrate(dnorm,-Inf,0)

0.5 with absolute error < 4.7e-05

integrate(dnorm,-Inf,100,abs.tol=0L)

0 with absolute error < 0

library(pracma)
integral(dnorm,-Inf,0)

[1] 0.5

integral(dnorm,-Inf,100,abstol=0)

[1] 0

这到底是怎么回事？我应该使用什么自适应方法？

Answer 1

查找 QAGI 和 QAGS 算法，似乎发生了以下情况:

1. 域 xϵ(-inf,b] 通过变换 x=b-(1-t)/t 从 t 进行映射，以便可以在 tϵ(0,1) 上计算积分。参见规范 here 。

2. adaptive quadrature algorithm被用来评估积分。将 limit=1 传递到您的 scipy 代码中会生成消息 “已达到最大分割数 (1)。” 传递 limit=2 不会产生此消息。这表明在算法的第 4 步中，积分的估计值 Q 和误差的估计值 ε 是相等的。

3. 这大概是因为估计中没有使用重要的点。在 tε(0,1] 的区间内使用 21 个均匀分布的点会产生范围为 80-100(大约)的 x 值。所有这些值都非常接近 0。算法中使用的值不是均匀分布的this page ，但大概会达到类似的结果。

因此，总而言之，从 (-inf,100] 到 (0,1] 的映射使积分估计中的采样值偏向 x=100 的终点。由于正态分布 pdf 是此处实际上为零，该算法不知道它缺少 x=0 附近的分布非零的区域，因此它永远不会分割以提高准确性。

此外，scipy 和 R 使用相同的算法，因此它们产生相同的结果是有道理的。

如果您从 -100 积分到 100，则中点 0 将是一个评估点，它允许算法按预期运行。但是，如果您从 -1000 积分到 100，算法会再次遗漏任何重要点，最终积分为 0。