scipy.integrate.quad对大数的精度

我尝试通过scipy.integrate.quad()计算这样的积分(实际上是pdf的指数分布的cdf)：

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中进行了讨论.

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

解决方法:

我相信问题是由于np.exp(-x)随着x的增加而迅速变小,由于数值精度有限而导致评估为零.例如,即使对于x小于x = 10 ** 2 *的x,np.exp(-x)的计算结果也为3.72007597602e-44,而10 ** 3或更高阶数的x值也为0.

我不知道Quad的实现细节,但是它可能会对给定的集成范围内要集成的功能执行某种采样.对于较大的积分上限,np.exp(-x)的大多数样本求值为零,因此积分值被低估了. (请注意,在这些情况下,所提供的绝对误差为四分之一,与积分值的阶次相同,这表明积分值不可靠.)

避免此问题的一种方法是将积分上限限制为一个数值,高于该数值数值函数会变得非常小(因此对积分值的贡献很小).从代码片段来看,值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]))

06002

我们还可以通过将精度提高到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]))

06004

通常,获得此精度的成本是增加的计算时间.

附注：毋庸置疑,如果您首先能够以分析方式评估您的积分(例如,在Sympy的帮助下),您会忘记上述所有内容.

标签：calculus,python,numpy,scipy
来源： https://codeday.me/bug/20191011/1895769.html