fork download
copy 
  1. # Incomplete gamma functions. (2.01)
  2. # (and chi-squared cumulative distribution)
  3. # @see Numerical Recipes in Fortran 77 (Ch. 6.2 Special Functions).
  4.  
  5. from math import lgamma, exp, log
  6.  
  7. # Regularized lower incomplete gamma function.
  8. # P(a, x) = γ(a, x) / Γ(a)
  9. def gammaP(a, x):
  10.     return 1 - gammaQ(a, x)
  11.  
  12. # Regularized upper incomplete gamma function.
  13. # Q(a, x) = Γ(a, x) / Γ(a)
  14. def gammaQ(a, x):
  15.     if x < 0 or a <= 0:
  16.         raise ValueError('bad arguments')
  17.     if x < a+1:
  18.         return 1 - p_series(a, x)
  19.     else:
  20.         return q_continued_fraction(a, x)
  21.  
  22. def p_series(a, x, itmax=1000, eps=3.0E-7):
  23.     if x == 0: # !
  24.         return 0.0
  25.     ap=a
  26.     sum_=1/a
  27.     del_=sum_
  28.     for _ in range(1, 1+itmax):
  29.         ap=ap+1
  30.         del_=del_*x/ap
  31.         sum_=sum_+del_
  32.         if abs(del_) < abs(sum_)*eps:
  33.             return sum_*exp(-x+a*log(x)-lgamma(a))
  34.     raise ValueError('a too large, itmax too small')
  35.  
  36. def q_continued_fraction(a, x, itmax=1000, eps=3.0E-7, fpmin=1.0E-30):
  37.     b=x+1-a
  38.     c=1/fpmin
  39.     d=1/b
  40.     h=d
  41.     for i in range(1, 1+itmax):
  42.         an=-i*(i-a)
  43.         b=b+2
  44.         d=an*d+b
  45.         if abs(d) < fpmin: d=fpmin
  46.         c=b+an/c
  47.         if abs(c) < fpmin: c=fpmin
  48.         d=1/d
  49.         del_=d*c
  50.         h=h*del_
  51.         if abs(del_-1) < eps:
  52.             return h*exp(-x+a*log(x)-lgamma(a))
  53.     raise ValueError('a too large, itmax too small')
  54.  
  55. # Chi-squared cumulative distribution function.
  56. def cdf(x, k):
  57.     return gammaP(k/2, x/2)
  58.  
  59. # Chi-squared p-value.
  60. def pvalue(x, k):
  61.     return 1 - cdf(x, k)
  62.  
  63. # Expect PV: [0.99, 0.95, 0.5, 0.05, 0.01]
  64. # @see Engineering Tables/Chi-squared Distribution - Wikibooks
  65. ks=[1, 2, 4, 10, 100, 1000]
  66. vs=[[1.57E-4, 0.00393,   0.455,    3.841,    6.635],
  67.     [ 0.0201,   0.103,   1.386,    5.991,    9.210],
  68.     [  0.297,   0.711,   3.357,    9.488,   13.277],
  69.     [  2.558,   3.940,   9.342,   18.307,   23.209],
  70.     [ 70.065,  77.929,  99.334,  124.342,  135.807],
  71.     [898.912, 927.594, 999.333, 1074.679, 1106.969]]
  72.  
  73. for k, v in zip(ks, vs):
  74.     print([round(pvalue(x, k), 2) for x in v])
Success #stdin #stdout 0.04s 10036KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[0.99, 0.95, 0.5, 0.05, 0.01]
[0.99, 0.95, 0.5, 0.05, 0.01]
[0.99, 0.95, 0.5, 0.05, 0.01]
[0.99, 0.95, 0.5, 0.05, 0.01]
[0.99, 0.95, 0.5, 0.05, 0.01]
[0.99, 0.95, 0.5, 0.05, 0.01]
https://ideone.com/vnAmvD
language:
Python 3 (python  3.9.5)
created:
5 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!