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 算法，似乎发生了以下情况：

域 xϵ(-inf,b] 通过变换 x=b-(1-t)/t 从 t 映射，以便可以在 tϵ(0, 1]. 查看规格 here.
正在使用adaptive quadrature algorithm计算积分。将 limit=1 传递到您的 scipy 代码中会产生消息 "The maximum number of subdivisions (1) has been achieved." 传递 limit=2 不会产生此消息。这表明在算法的第 4 步中，积分的估计值 Q 和误差的估计值 ε 相等。
这大概是因为估计中没有使用重要的点。在 tε(0,1] 的区间内使用 21 个均匀分布的点得到 x-values，范围从 80-100（大约）。所有这些值都非常接近 0。算法中使用的值不均匀每 this page 间隔一次，但大概会达到类似的结果。

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

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

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

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

The functions "integrate.quad" in python and "integral",'integrate' in R give erroneous results

python

r

integrate