利用积分法推导CDF的精度

Question

我有两种方法来推导一个正态分布的随机变量在一个区间内的概率。第一个也是最直接的是:

import scipy.stats
print scipy.stats.norm.cdf(6) - scipy.stats.norm.cdf(5)
# 2.85664984223e-07

第二种是通过整合pdf:

import scipy.integrate
print scipy.integrate.quad(scipy.stats.norm.pdf, 5, 6)[0]
# 2.85664984234e-07

这种情况下的差异非常小，但这并不意味着它不会在其他分布或积分限制下变得更大。你能说出哪个更准确吗?为什么?

顺便说一下，第一种方法似乎至少快10倍，所以如果它也更准确(这是我的猜测，因为它有点专业化)，那么它就是完美的。

Answer 1

在这种特殊情况下，给定这些特定的数字，quad方法实际上会更准确。当然，CDF本身可以快速准确地计算出来，但是看看实际的数字:

>>> scipy.stats.norm.cdf(6), scipy.stats.norm.cdf(5)
(0.9999999990134123, 0.99999971334842808)

当你差两个非常相似的量时，你会失去准确性。在集成过程中，如果编码员小心处理他们的求和，类似的问题可以得到一定程度的缓解。

无论如何，我们可以使用mpmath对高分辨率计算进行检查:

>>> via_cdf = scipy.stats.norm.cdf(6)-scipy.stats.norm.cdf(5)
>>> via_quad = scipy.integrate.quad(scipy.stats.norm.pdf, 5, 6)[0]
>>> import mpmath
>>> mpmath.mp.dps = 100
>>> def cdf(x): return 0.5 * (1 + mpmath.erf(x/mpmath.sqrt(2)))
>>> highres = cdf(6)-cdf(5)
>>> highres
mpf('0.0000002856649842341562135330514687422473118357532223619105443630157837185833042478210791954518847897468442097')
>>> float((highres - via_quad)/highres)
-2.3824773334590333e-16
>>> float((highres - via_cdf)/highres)
3.86659439572868e-11

Answer 2

第一个调用scipy.special中包含的cdf的实现。后者实际上做了积分。前者可能更准确(因为它只受限于计算机评估CDF的能力，而不受数值积分引入的任何误差的限制)。在实践中，除非您需要的结果好于小数点后6位，否则您可能会很好。