使用scipy.integrate.quad时结果不连续 [英] Discontinuity in results when using scipy.integrate.quad
问题描述
使用scipy.integrate.quad时，我发现了一个奇怪的行为.此行为也出现在Octave的quad函数中，这使我相信它可能与QUADPACK本身有关.有趣的是，使用完全相同的Octave代码，此行为在MATLAB中不会出现.
I've discovered a strange behavior when using scipy.integrate.quad. This behavior also shows up in Octave's quad function, which leads me to believe that it may have something to do with QUADPACK itself. Interestingly enough, using the exact same Octave code, this behavior does not show up in MATLAB.
关于这个问题.我在数值上整合了各种边界上的对数正态分布.因为F是对数正态的cdf，a是下限，b是上限，所以我发现在某些情况下，
On to the question. I'm numerically integrating a lognormal distribution over various bounds. For F is cdf of lognormal, a is lower bound and b is upper bound, I find that under some conditions,
当b是一个非常大的数字"时，
integral(F，a，b)= 0，而
integral(F, a, b) = 0 when b is a "very large number," while
integral(F，a，b)=当b为np.inf时的正确极限. (或者只是Inf for Octave.)
integral(F, a, b) = the correct limit when b is np.inf. (or just Inf for Octave.)
下面是一些示例代码，用于演示操作:
Here's some example code to show it in action:
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)评估为近似正确的值，但是对于所有值832< = b<，quad(F，1.0，b)变为零. np.inf?
What is going on here? Why does quad(F, 1.0, 500) evaluate to approximately the correct thing, but quad(F, 1.0, b) goes to zero for all values 832 <= b < np.inf?
推荐答案
虽然我对QUADPACK并不十分熟悉，但是自适应集成通常通过提高分辨率来起作用，直到答案不再改善为止.对于大多数时间间隔(使用F(10)==9.356e-116)，您的函数是如此接近0，以至于Quad选择的初始网格点的改进可以忽略不计，并且它决定积分必须接近0.基本上，如果数据隐藏在积分范围的很小的子区间中，quad最终将找不到它.
While I'm not exactly familiar with QUADPACK, adaptive integration generally works by increasing resolution until the answer no longer improves. Your function is so close to 0 for most of the interval (with F(10)==9.356e-116) that the improvement is negligible for the initial grid points that quad chooses, and it decides that the integral must be close to 0. Basically, if your data hides in a very narrow subinterval of the range of integration, quad eventually won't be able to find it.
对于从0到inf的积分，显然不能将间隔细分为有限数量的间隔，因此quad在计算积分之前需要进行一些预处理.例如，变量y=1/(1+x)的更改会将区间0..inf映射到0..1.将该间隔细分会从原始函数中采样更多接近零的点，从而使quad能够找到您的数据.
For integration from 0 to inf, the interval obviously cannot be subdivided into a finite number of intervals, so quad will need some preprocessing before computing the integral. For example, a change of variables like y=1/(1+x) would map the interval 0..inf to 0..1. Subdividing that interval will sample more points near zero from the original function, enabling quad to find your data.
这篇关于使用scipy.integrate.quad时结果不连续的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
