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。