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，绝对误差< 4.7e-05

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

0，绝对误差 < 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] 上计算积分。请参阅此处的规格 .

2. 自适应求积算法用于评估积分。将

   limit=1

   传递到您的

   scipy

   代码中会产生消息 “已达到最大细分数 (1)。” 传递

   limit=2

   不会产生此消息。这表明在算法的第 4 步中，积分 Q 的估计值和误差 ε 的估计值是相等的。

3. 发生这种情况可能是因为估计中没有使用重要点。在 tϵ(0,1] 区间内使用 21 个均匀间隔的点会产生范围为 80-100（大约）的 x 值。所有这些值都非常接近 0。算法中使用的值不是均匀间隔的此页面，但大概可以达到类似的结果。

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

此外，

scipy

和

R

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

如果从-100到100积分，中点0将成为评估点，这使得算法能够按预期运行。但是，如果您从 -1000 积分到 100，算法会再次错过任何重要点，最终积分为 0。