使用Scipy计算高斯函数尾部积分时结果为零而非8.19e-26

Question

我正在尝试对一个高斯函数进行积分，但积分的范围在高斯函数的尾部，所以用integrate.quad计算出来的结果是零。有没有办法对这个高斯函数进行积分，得到一个非常小的结果呢？

这个函数的被积函数是：

sigma = 9.5e-5
integrand = lambda delta: (1./(np.sqrt(2*np.pi)*sigma))*np.exp(-(delta**2)/(2*sigma**2))

我需要在 10^-3 到 0.3 之间进行积分。

通过Wolfram Alpha，我得到了一个结果是8.19e-26，但用Scipy的Romberg积分方法却得到了零。有没有办法调整Scipy的参数，以便对这样的小结果进行积分呢？

Answer 1

谢谢你的帮助，

经过更多的讨论，我选择了数值积分这个方法。在用C++脚本检查后，我发现如果我在scipy.integrate.romberg中设置divmax = 120，就能得到和Wolfram Alpha一样的结果。

不过，这个方法计算起来花费的时间很长。我会尝试使用误差函数，看看能否理解它。

祝好

Answer 2

是的。

>>> from scipy.special import erfc
>>> erfc(1e-3/9.5e-5/np.sqrt(2.))
6.534205277034387e-26

在这个尾部的部分，使用补充误差函数（erfc）会更好，或者你也可以考虑使用erfcx，它是补充误差函数经过exp(x**2)缩放后的版本。

Answer 3

让我们先定义一下，F(x; s) 是一个正态分布（也叫高斯分布）的累积分布函数，标准差是 s。你现在要计算的是 F(x1;s) - F(x0;s)，其中 x0 = 1e-3 和 x1 = 0.3。

这个计算可以改写成 S(x0;s) - S(x1;s)，这里的 S(x;s) = 1 - F(x;s) 被称为“生存函数”。

你可以使用 scipy.stats 中的 norm 对象的 sf 方法来计算这个值。

In [99]: x0 = 1e-3

In [100]: x1 = 0.3

In [101]: s = 9.5e-5

In [102]: from scipy.stats import norm

In [103]: norm.sf(x0, scale=s)
Out[103]: 3.2671026385171459e-26

In [104]: norm.sf(x1, scale=s)
Out[104]: 0.0

需要注意的是，norm.sf(x1, scale=s) 的结果是 0。这个表达式的确切值是一个非常小的数字，甚至小到无法用 64 位浮点数表示（正如 @Zhenya 在评论中提到的）。

所以这个计算的结果是 3.267e-26。

你也可以使用 scipy.special.ndtr 来计算这个值。ndtr 计算的是标准正态分布的累积分布函数，并且由于对称性，S(x; s) = ndtr(-x/s)。

In [105]: from scipy.special import ndtr

In [106]: ndtr(-x0/s)
Out[106]: 3.2671026385171459e-26

---

如果你想通过数值积分来得到相同的结果，你需要尝试调整积分算法的误差控制参数。例如，使用 scipy.integrate.romberg 时，我调整了 divmax 和 tol，如下所示：

In [60]: from scipy.integrate import romberg

In [61]: def integrand(x, s):
   ....:     return np.exp(-0.5*(x/s)**2)/(np.sqrt(2*np.pi)*s)
   ....: 

In [62]: romberg(integrand, 0.001, 0.3, args=(9.5e-5,), divmax=20, tol=1e-30)
Out[62]: 3.2671026554875259e-26

使用 scipy.integrate.quad 时，我需要告诉它 0.002 是一个“特殊”的点，这样它才能多花点力气处理：

In [81]: from scipy.integrate import quad

In [82]: p, err = quad(integrand, 0.001, 0.3, args=(9.5e-5,), epsabs=1e-32, points=[0.002])

In [83]: p
Out[83]: 3.267102638517144e-26

In [84]: err
Out[84]: 4.769436484142494e-37