在很长的时间间隔内积分高斯答案

【问题标题】：Integrating a gaussian over a very long interval在很长的时间间隔内积分高斯
【发布时间】：2019-05-15 09:27:08
【问题描述】：

我想在一个非常大的区间内积分一个高斯函数。我选择了spicy.integrate.quad 函数进行集成。该功能似乎仅在我选择足够小的间隔时才起作用。当我使用下面的代码时，

from scipy.integrate import quad
from math import pi, exp, sqrt

def func(x, mean, sigma):
    return 1/(sqrt(2*pi)*sigma) * exp(-1/2*((x-mean)/sigma)**2) 

print(quad(func, 0, 1e+31, args=(1e+29, 1e+28))[0]) # case 1
print(quad(func, 0, 1e+32, args=(1e+29, 1e+28))[0]) # case 2
print(quad(func, 0, 1e+33, args=(1e+29, 1e+28))[0]) # case 3
print(quad(func, 1e+25, 1e+33, args=(1e+29, 1e+28))[0]) # case 4

然后打印以下内容。

1.0
1.0000000000000004
0.0
0.0

为了获得合理的结果，我不得不尝试多次更改积分的下限/上限，并凭经验将其确定为 [0, 1e+32]。这对我来说似乎很冒险，因为当高斯函数的均值和西格玛发生变化时，我总是不得不尝试不同的界限。

有没有一种明确的方法可以将函数从 0 集成到 1e+50 而不用担心边界？如果不是，您从一开始就期望哪些边界会给出非零值？

【问题讨论】：

查看我对这个问题的回答：stackoverflow.com/questions/29179778/…
也相关：stackoverflow.com/questions/47193045/…

标签： python scipy integral

【解决方案1】：

简而言之，你不能。

在这个长间隔上，高斯非零的区域很小，在integrate.quad 引擎下工作的自适应过程看不到它。除非偶然，几乎任何自适应程序都会如此。

【讨论】：

【解决方案2】：

注意，

正常随机变量的 CDF 称为ϕ(x)，因为它不能用elementary function 表示。所以采取ϕ((b-m)/s) - ϕ((a-m)/s)。另请注意ϕ(x) = 1/2(1 + erf(x/sqrt(2)))，因此您无需调用.quad 来实际执行集成，并且使用scipy 中的erf 可能会有更好的运气。

from scipy.special import erf

def prob(mu, sigma, a, b):
    phi = lambda x: 1/2*(1 + erf((x - mu)/(sigma*np.sqrt(2))))
    return phi(b) - phi(a)

这可能会给出更准确的结果（比上面的）

>>> print(prob(0, 1e+31, 0, 1e+50))
0.5
>>> print(prob(0, 1e+32, 1e+28, 1e+29))
0.000359047985937333
>>> print(prob(0, 1e+33, 1e+28, 1e+29))
3.5904805169684195e-05
>>> print(prob(1e+25, 1e+33, 1e+28, 1e+29))
3.590480516979522e-05

并避免您遇到的强烈的floating point 错误。但是，您整合的区域太小了，您仍然可以看到0。

【讨论】：