将高斯函数的尾部与 Scipy 集成，给出零而不是 8.19e-26答案

【问题标题】：Integration of the tail of a Gaussian function with Scipy, giving zero instead of 8.19e-26将高斯函数的尾部与 Scipy 集成，给出零而不是 8.19e-26
【发布时间】：2014-03-25 10:03:43
【问题描述】：

我正在尝试对高斯函数进行积分，限制在高斯尾部内部，所以尝试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 中的旋钮来整合这么小的结果吗？

【问题讨论】：

尝试 romberg 与参数 tol=1e-28, rtol=1e-28;你会得到3.25e-25这与mathematica结果一致
您能否详细说明您是如何从 Wolfram Alpha 获得 8.19e-26 的？

标签： python python-2.7 numpy scipy

【解决方案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

【讨论】：

正如我对@Zhenya 回答的评论所说，这与 Wolfram Alpha 的结果相去甚远。尝试在 Wolfram Alpha 中评估 SurvivalFunction[NormalDistribution[0, 9.5*10^(-5)], .003] 或 CDF[NormalDistribution[0, 9.5*10^(-5)], -.003]。结果是O(e-219)。
但是下限是0.001，不是0.003。
当然，我的错。这样就解决了。

【解决方案2】：

是的。

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

你最好使用补足误差函数 (erfc) 或者erfcx，它是由exp(x**2) 缩放的补足误差函数。

【讨论】：

您正在测试集成的“简单”极端，即delta=1e-3。如果您现在为delta=.3 运行测试，您的测试将失败，即erfc(.3/9.5e-5/np.sqrt(2.))。我不会对你投反对票，因为错误很容易发生。
在上面的评论中是erfc(.3/9.5e-5/np.sqrt(2.))==0，我错过了==0 位
@flebool 鉴于erfcx(0.3/9.5e-5 / np.sqrt(2.)) 大致为2e-4，您说的是1e-2165457 顺序的附加校正
我同意，这实在是太小了。但即使是另一个极端也非常小，在数值噪声范围内。 Wolfram Alpha 的结果是10^(-219) 的顺序，见link
我这里的意思是上面的结果erfc(1e-3/9.5e-5/np.sqrt(2.))=6.534205277034387e-26是错误的，如Wolfram Alpha的计算所示。

【解决方案3】：

感谢您的帮助，

经过更多咨询后，我转到数值积分选项，在检查 c++ 脚本后，我发现如果我在 scipy.integrate.romberg 设置 divmax = 120，我会得到与 Wolfram Alpha 相同的结果.

但是这个解决方案需要大量时间来计算。我将尝试使用错误函数，看看我是否能理解它..

干杯

【讨论】：

这是一个标准的、定义明确的、广泛使用的数量。您不应该使用数值积分来计算它，特别是因为数值积分给出了错误的答案。 @flebool 给出了这个值的 Mathematica 公式：SurvivalFunction[NormalDistribution[0, 9.5*10^(-5)], .001]。在 Wolfram Alpha 中尝试一下。如果您仍然不相信，请尝试通过将两个子区间（从 0.001 到 0.002 以及从 0.002 到 0.3）上的数值积分结果相加来计算数值积分。
所以我尝试了 norm.sf 解决方案，它确实给出了相同的答案，而且速度更快。