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

Question

使用scipy.integrate.quad时，我发现了一个奇怪的行为。此行为也出现在Octave的quad函数中，这使我相信它可能与QUADPACK本身有关。有趣的是，使用完全相同的Octave代码，此行为未在MATLAB中显示。

关于这个问题。我在数值上整合了各种边界上的对数正态分布。因为F是对数正态的cdf，a是下限，b是上限，所以我发现在某些情况下，

当b是一个“非常大的数字”时，integral（F，a，b）= 0，而

当b为np.inf时，integrate（F，a，b）=正确的极限。（或者只是Inf for Octave。）

以下是一些示例代码，以展示其实际效果：

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 

这里发生了什么？为什么对所有值832 <= b <np.inf，quad（F，1.0，500）求近似近似正确的值，但是quad（F，1.0，b）变为零？

Answer 1

尽管我对QUADPACK并不十分熟悉，但是自适应集成通常通过提高分辨率直到效果不再改善而起作用。您的函数在大部分时间间隔内（带有F(10)==9.356e-116）都非常接近0，以至于Quad选择的初始网格点的改善可以忽略不计，并且它决定积分必须接近0。基本上，如果数据隐藏在集成范围的较小子间隔，quad最终将无法找到它。

对于从0到的积分inf，显然不能将间隔细分为有限数量的间隔，因此quad在计算积分之前需要进行一些预处理。例如，如变量的更改y=1/(1+x)会将间隔映射0..inf到0..1。细分该间隔将从原始函数中采样更多接近零的点，从而quad找到您的数据。