大数上的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的使用解决了这个问题，但所有使用上限集成限制的尝试都得到了不正确的答案。

在#5428在GitHub上发行中讨论了类似的情况。

在集成其他密度函数时，我应该如何避免这样的错误？

Answer 1

我相信这个问题是由于np.exp(-x)随着x的增加而迅速变得很小，这导致了由于有限的数值精度而被评定为零。例如，即使对于像x这样小的x=10**2*，np.exp(-x)的计算结果也是3.72007597602e-44，而order 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])