scipy.integrate.quad 大数精度

Question

我尝试通过以下方式计算这样的积分（实际上是指数分布的 cdf 及其 pdf）scipy.integrate.quad()：

import numpy as np
from scipy.integrate import quad

def g(x):
    return .5 * np.exp(-.5 * x)

print quad(g, a=0., b=np.inf)
print quad(g, a=0., b=10**6)
print quad(g, a=0., b=10**5)
print quad(g, a=0., b=10**4)

结果如下：

(1.0, 3.5807346295637055e-11)
(0.0, 0.0)
(3.881683817604194e-22, 7.717972744764185e-22)
(1.0, 1.6059202674761255e-14)

尽管使用np.inf解决了问题，但所有尝试使用大的积分上限都会产生错误的答案。

GitHub 上的 scipy issue #5428 中讨论了类似的情况。

我应该怎么做才能避免在集成其他密度函数时出现这种错误？

Answer 1

我相信这个问题是由于np.exp(-x)随着x增加而迅速变得非常小，这导致由于数值精度有限而评估为零。例如，即使x小到x=10**2*，也np.exp(-x)计算为3.72007597602e-44，而xorder10**3或更高的值会导致0。

我不知道 的实现细节quad，但它可能会在给定的集成范围内对要集成的函数执行某种采样。对于较大的积分上限，大多数样本np.exp(-x)评估为零，因此低估了积分值。（请注意，在这些情况下，提供的绝对误差 byquad与积分值具有相同的顺序，这是后者不可靠的指标。）

避免此问题的一种方法是将积分上限限制为一个值，高于该值，数值函数变得非常小（因此，对积分值的贡献很小）。从您的代码片段来看， 的值10**4似乎是一个不错的选择，但是， 的值10**2也会导致对积分的准确评估。

避免数值精度问题的另一种方法是使用以任意精度算术执行计算的模块，例如mpmath. 例如， for x=10**5，mpmath计算exp(-x)如下（使用原生mpmath指数函数）

import mpmath as mp
print(mp.exp(-10**5))

3.56294956530937e-43430

注意这个值有多小。使用标准硬件数值精度（由 使用numpy），该值变为0.

mpmath提供了一个积分函数 ( mp.quad)，它可以为积分上限的任意值提供积分的准确估计。

import mpmath as mp

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))

1.0
0.999999650469474
0.999999999996516
0.999999999999997

我们还可以通过将精度提高到50小数点（15标准精度）来获得更准确的估计

mp.mp.dps = 50; 

print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, mp.inf]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**13]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**8]))
print(mp.quad(lambda x : .5 * mp.exp(-.5 * x), [0, 10**5]))

1.0
0.99999999999999999999999999999999999999999829880262
0.99999999999999999999999999999999999999999999997463
0.99999999999999999999999999999999999999999999999998

一般来说，获得这种精度的代价是增加了计算时间。

PS：不用说，如果您能够首先分析地评估您的积分（例如，在 的帮助下Sympy），您可以忘记以上所有内容。