使用 scipy.integrate.quad 时，正函数的增加界限会减少积分！怎么解决？

Question

我正在使用 scipy.integrate 来计算定积分。我编写了以下代码片段来评估概率密度函数 (PDF) 和累积分布函数 (CDF)，以评估 PDF 中的 CDF，如下所示：

inf = 10000

def f_normal(m, tau): #normal pdf with mean np.sqrt(m- 1/2) and var 1/pi
    tmp = 1/np.sqrt(np.pi) * np.exp(  -  ( tau - np.sqrt(m- 1/2) )**2   )
    return tmp

def F_normal(m, x): #integral of f_normal from -inf to x
    return integrate.quad(lambda tau: f_normal(m, tau), -inf, x)

但是当我评估时：

m=10    
print(f_normal(m, 10))    
print(F_normal(m, inf))

我得到

(0.0, 0.0)

换

F_normal(m, inf)

，而我真的很期待1。我不明白发生了什么事。我还尝试理解命令integrate.quad，因此我得到了这些奇怪的值：虽然增加界限会减少正函数PDF的定积分，因此没有意义：

integrate.quad(lambda tau: f_normal(m, tau), -10, 10)
Out[30]: (1.0, 2.2977017217085778e-09)

integrate.quad(lambda tau: f_normal(m, tau), -100, 100)
Out[31]: (1.0, 8.447041410391393e-09)

integrate.quad(lambda tau: f_normal(m, tau), -1000, 1000)
Out[32]: (0.9999934640807174, 1.857441875165816e-09)

integrate.quad(lambda tau: f_normal(m, tau), -10000, 10000)
Out[33]: (5.659684283308144e-150, 1.1253180602819657e-149)

integrate.quad(lambda tau: f_normal(m, tau), -10000, 100)
Out[34]: (0.0, 0.0)

integrate.quad(lambda tau: f_normal(m, tau), -10000, 1000000)
Out[35]: (0.0, 0.0)

我有点困惑：感谢帮助！

Answer 1

问题是

inf

太大了。您正在尝试集成此功能，

但是您要发送到

integrate.quad

的是，

求积例程不太可能在值本质上不为 0 的任何地方对该函数进行采样，从而导致错误的结果。

幸运的是，使用

np.inf

，您可以告诉算法积分在无限区间内，并且它会更智能地查找函数的相关行为在哪里。此代码给出了正确的积分：

def f_normal(m, tau): #normal pdf with mean np.sqrt(m- 1/2) and var 1/pi
    tmp = 1/np.sqrt(np.pi) * np.exp(  -  ( tau - np.sqrt(m- 1/2) )**2   )
    return tmp

def F_normal(m, x): #integral of f_normal from -inf to x
    return integrate.quad(lambda tau: f_normal(m, tau), -np.inf, x)

m=10    
print(f_normal(m, m-.5))    
print(F_normal(m, np.inf))

或者，您可以在函数不太接近 0 的地方对较小的有限区间进行积分。在这种情况下，中心 (

np.sqrt(m-1/2)

) 正负 100 是可靠的（其中长度为 20,000 的原始区间太大）：

def f_normal(m, tau): #normal pdf with mean np.sqrt(m- 1/2) and var 1/pi
    tmp = 1/np.sqrt(np.pi) * np.exp(  -  ( tau - np.sqrt(m- 1/2) )**2   )
    return tmp

def F_normal(m, x): #integral of f_normal from -inf to x
    return integrate.quad(lambda tau: f_normal(m, tau), np.sqrt(m- 1/2)-100, x)

m=10000
print(f_normal(m, m-.5))    
print(F_normal(m, np.sqrt(m- 1/2)+100))