cdf精度scipy.stats公司

Question

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

对于给定的区间，我需要将这个分布估计为一个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函数精确？在

Answer 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），走

Answer 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）：

^{pr2}$

使用quad来获得任何需要的规范化应该很简单。在