python - 使用 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)

我有点困惑：感谢帮助！

标签： pythonscipyintegralintegrate

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))