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