使用scipy.integrate.quad时结果的不连续性

Question

我发现使用scipy.integrate.quad的时候有个奇怪的现象。这种现象在Octave的quad函数中也出现过，这让我觉得可能和QUADPACK本身有关。有趣的是，使用完全相同的Octave代码，这种现象在MATLAB中却没有出现。

接下来说说我的问题。我正在对一个对数正态分布进行数值积分，积分的上下限分别是a和b，F是对数正态分布的累积分布函数。在某些情况下，我发现：

当b是一个“非常大的数字”时，

integral(F, a, b) = 0；而当b是np.inf（或者在Octave中就是Inf）时，

integral(F, a, b) = 正确的极限。

这里有一些示例代码来展示这个现象：

from __future__ import division
import numpy as np
import scipy.stats as stats
from scipy.integrate import quad

# Set up the probability space:
sigma = 0.1
mu = -0.5*(sigma**2)  # To get E[X] = 1
N = 7
z = stats.lognormal(sigma, 0, np.exp(mu))


# Set up F for integration:
F = lambda x: x*z.pdf(x)

# An example that appears to work correctly:
a, b = 1.0, 10
quad(F, a, b)
# (0.5199388..., 5.0097567e-11)

# But if we push it higher, we get a value which drops to 0:
quad(F, 1.0, 1000)
# (1.54400e-11, 3.0699e-11)

# HOWEVER, if we shove np.inf in there, we get correct answer again: 
quad(F, 1.0, np.inf)
# (0.5199388..., 3.00668e-09)


# If we play around we can see where it "breaks:"
quad(F, 1.0, 500)   # Ok
quad(F, 1.0, 831)   # Ok
quad(F, 1.0, 832)   # Here we suddenly hit close to zero. 
quad(F, 1.0, np.inf)   # Ok again 

这到底是怎么回事呢？为什么quad(F, 1.0, 500)的结果大致是正确的，而quad(F, 1.0, b)在所有832 <= b < np.inf的情况下都变成了零呢？

Answer 1

试着降低错误容忍度

>>> quad(F, a, 1000, epsabs=1.49e-11)
(0.5199388058383727, 2.6133800952484582e-11)

我觉得数值积分对某些设置特别敏感。你可以通过调用 quad(..., full_output=1) 来调试它，然后仔细分析详细输出的内容。如果这个回答不够满意，我很抱歉。

Answer 2

虽然我对QUADPACK不太熟悉，但自适应积分的工作原理一般是通过不断提高精度，直到结果不再有明显改善。你的函数在大部分区间内都非常接近0（比如F(10)==9.356e-116），所以对于quad选择的初始网格点来说，改进几乎可以忽略不计，因此它认为这个积分应该也接近0。简单来说，如果你的数据藏在积分范围内一个非常狭窄的子区间里，quad最终就找不到它了。

对于从0到inf的积分，这个区间显然不能被划分成有限个小区间，所以在计算积分之前，quad需要一些预处理。比如，像y=1/(1+x)这样的变量变换，可以把区间0..inf映射到0..1。这样划分后，能从原始函数中在接近0的地方取到更多的点，这样quad就能找到你的数据了。