python - 将高斯函数的尾部与 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

设 F(x; s) 为正态(即高斯)分布的 CDF，其中标准差 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 计算标准正态分布的 CDF，根据对称性，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