scipy.integrate.quad在处理大数时的精度问题

Question

scipy.integrate.quad在处理大数时的精度问题

4

我试图通过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问题＃5428中讨论。

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

- Denis Korzhenkov

2个回答

3

使用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])

- pv.

与 np.quad(g, a=0., b=100) 的输出相比较，这种方法似乎将上限设置为100，而不考虑实际用户输入。当然，对于 OP 的目的来说，这可能是完全可以接受的。 - Stelios

它并不会。积分器确实会对 x > 100 的函数进行采样，但这当然是一个基本事实，即该积分的那部分只会给出非常小的贡献。 - pv.

@pv，在阅读quad文档字符串后，我不明白你的建议如何有帮助。1和100点不是不连续点。 - Denis Korzhenkov

1

@DenisKorzhenkov：它强制积分器在这些点上采样函数。否则，对于大的积分区间，它将在a + eps * (b - a)这些点采样 --- 但如果b-a非常大，它将错过接近x=0的峰值。 - pv.

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Stelios · Accepted Answer

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

我不知道quad的具体实现，但它可能执行一些函数的采样以在给定的积分范围内进行积分。对于较大的上限积分，np.exp(-x)的大多数样本评估为零，因此低估了积分值。（请注意，在这些情况下，quad提供的绝对误差与积分值的数量级相同，这表明后者不可靠。）

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

避免数值精度问题的另一种方法是使用一个执行任意精度算术计算的模块，例如mpmath。例如，对于x=10**5，mpmath使用本地mpmath指数函数评估exp(-x)如下：

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）评估积分，那么您可以忘记上述所有内容。