通常每当我遇到精度问题时,我使用的第一个工具就是 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 来获得所需的任何归一化。