scipy.integrate.quad 大数精度

Question

我尝试通过 scipy.integrate.quad():

计算这样一个积分（实际上是指数分布的 cdf 及其 pdf）

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 解决了问题，但所有使用大积分上限的尝试都会产生不正确的答案。

scipy issue #5428 at GitHub 中讨论了类似的案例。

如何避免在积分其他密度函数时出现这样的错误？

Answer 1

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

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

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

另一种避免数值精度问题的方法是使用以任意精度算法执行计算的模块，例如 mpmath。例如，对于 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

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

P.S.: 不言而喻，如果你能够首先分析地评估你的积分（例如，在 Sympy 的帮助下），你可以忘记以上所有内容。

Answer 2

使用 points 参数告诉算法你的函数的支持大致是：

import numpy as np
from scipy.integrate import quad

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

print quad(g, a=0., b=10**3, points=[1, 100])
print quad(g, a=0., b=10**6, points=[1, 100])
print quad(g, a=0., b=10**9, points=[1, 100])
print quad(g, a=0., b=10**12, points=[1, 100])

scipy.integrate.quad 大数精度

scipy.integrate.quad precision on big numbers

python

numpy

scipy

calculus