通过积分推导CDF的准确性

Question

我有两种方法来计算一个正态分布的随机变量落在某个区间内的概率。第一种方法是最简单直接的：

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

第二种方法是通过对概率密度函数进行积分：

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

在这种情况下，两者的差别非常小，但这并不意味着在其他分布或积分范围内差别不会变大。你能告诉我哪种方法更准确，以及为什么吗？

顺便提一下，第一种方法似乎至少快了10倍，所以如果它也更准确（我猜是这样，因为它有点专业），那就太完美了。

Answer 1

第一个调用的是在 scipy.special 中包含的累积分布函数（cdf）的实现。第二个实际上是进行积分计算。第一个方法可能更准确，因为它只受限于计算机评估累积分布函数的能力，而不是由于数值积分引入的错误。在实际应用中，除非你需要结果精确到小数点后六位，否则你用第一个方法就足够了。

Answer 2

在这个特定的情况下，考虑到这些特定的数字，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