在很长的时间间隔内对高斯进行积分

Question

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

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

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

是否有一种清晰的方法可以在不考虑边界的情况下将函数从 0 集成到 1e+50？如果不是，您如何从一开始就期望哪些边界会给出非零值？

Answer 1

简而言之，你不能。

在这个长间隔上，高斯 non-zero 的区域很小，在 integrate.quad 引擎盖下工作的自适应程序无法看到它。除非偶然，否则几乎任何适应性程序都会如此。

Answer 2

通知，

正常随机变量的 CDF 称为 ϕ(x)，因为它不能用 scipy 中的 elementary function. So take ϕ((b-m)/s) - ϕ((a-m)/s). Also note that ϕ(x) = 1/2(1 + erf(x/sqrt(2))) so you need not call .quad to actually perform an integration and may have better luck with 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.

在很长的时间间隔内对高斯进行积分

Integrating a gaussian over a very long interval

python

scipy

integral