scipy.stats 中 cdf 的精度答案

【问题标题】：Precision of cdf in scipy.statsscipy.stats 中 cdf 的精度
【发布时间】：2011-06-09 19:23:43
【问题描述】：

我使用chi2 分布作为模拟系统的理论问题。

对于给定的区间，我需要将此分布估计为 PMF，该 PMF 定义为该区间内 PDF 的积分。该值应接近区间中心的 PDF 值，但可能略有不同，具体取决于 PDF 的形状。

这是我的工作：

import numpy
from scipy.stats import chi2

dist = chi2(10)
nbins = 120

F = dist.cdf(numpy.arange(nbins+1))
pmf = F[1:] - F[:-1] # surface inside the interval
pmf /= pmf.sum() # Normalisation

问题是chi2.cdf(100, 10) 及以上给出的正好是 1.0。所以我能得到的最小值是 1.11e-16 左右。但是chi2.pdf(100, 10) 不完全是 0（大约是 2.5e-17）。

我的问题是：我怎样才能更精确地获得我的 pmf 估计（可能高达 1e-25）？为什么 cdf 函数的精度不如 pdf 函数？

【问题讨论】：

标签： python numpy scipy precision cdf

【解决方案1】：

cdf 在等于 1 的浮点精度范围内，但 sf 接近于零，因此微小的差异 1e-20 不会被大 1 覆盖。（请参阅 JABS 参考）

>>> probs_from_cdf = np.diff(stats.chi2.cdf(np.arange(nbins+1), 10))
>>> probs_from_sf = np.diff(stats.chi2.sf(np.arange(nbins+1)[::-1], 10))[::-1]
>>> probs_from_sf[:4]
array([ 0.00017212,  0.00348773,  0.01491609,  0.03407708])
>>> probs_from_cdf[:4]
array([ 0.00017212,  0.00348773,  0.01491609,  0.03407708])
>>> probs_from_cdf[-5:]
array([ 0.,  0.,  0.,  0.,  0.])
>>> probs_from_sf[-5:]
array([  1.94252577e-20,   1.21955220e-20,   7.65430774e-21,
         4.80270079e-21,   3.01259913e-21])

我不知道 sf 的准确范围，即 scipy.special.chdtrc(df, x) 有多远

【讨论】：

谢谢！这正是我一直在寻找的。作为奖励，你向我展示了我不知道的 diff 函数。

【解决方案2】：

通常，每当我遇到精度问题时，我使用的第一个工具是 mpmath。 90% 的时间它只是工作（tm），足够快。在这种情况下，我们可以这样写：

import mpmath
mpmath.mp.dps = 50 # decimal digits of precision

def pdf(x,k):
    x,k = mpmath.mpf(x), mpmath.mpf(k)
    if x < 0: return 0
    return 1/(2**(k/2) * mpmath.gamma(k/2)) * (x**(k/2-1)) * mpmath.exp(-x/2)

def cdf(x,k): 
    x,k = mpmath.mpf(x), mpmath.mpf(k) 
    return mpmath.gammainc(k/2, 0, x/2, regularized=True)

def cdf_via_quad(s,k):
    return mpmath.quad(lambda x: pdf(x,k), [0, s])

给予（使用你的 F）：

>>> pdf(2,10)
mpf('0.0076641550244050483665734118783637680717877318964951605')
>>> cdf(2,10)
mpf('0.003659846827343712345456455812710150667594853455628779')
>>> cdf_via_quad(2,10)
mpf('0.003659846827343712345456455812710150667594853455628779')
>>> F[2]
0.0036598468273437131
>>> pdf(100,10)
mpf('2.5113930312030179466371651256862142900427508479560716e-17')
>>> cdf(100,10)
mpf('0.99999999999999994550298017079470664906667698474760744')
>>> cdf_via_quad(100,10)
mpf('0.99999999999999994550298017079470664906667698474760744')
>>> F[100]
1.0

应该可以直接使用 quad 来获得所需的任何归一化。

【讨论】：