1
2
3
4
5
6
7
8
9
10 from math import sqrt, log, exp
11
13 if df < 1:
14 raise ValueError, "df must be at least 1"
15 if stat < 0:
16 raise ValueError, "The test statistic must be positive"
17 x = 0.5 * stat
18 alpha = df / 2.0
19 prob = 1 - _incomplete_gamma(x, alpha)
20 return prob
21
23 """Compute the log of the gamma function for a given alpha.
24
25 Comments from Z. Yang:
26 Returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.
27 Stirling's formula is used for the central polynomial part of the procedure.
28 Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
29 Communications of the Association for Computing Machinery, 9:684
30 """
31 if alpha <= 0:
32 raise ValueError
33 x = alpha
34 f = 0
35 if x < 7:
36 f = 1
37 z = x
38 while z<7:
39 f *= z
40 z += 1
41 x = z
42 f = -log(f)
43 z = 1 / (x * x)
44 return f + (x-0.5)*log(x) - x + .918938533204673 \
45 + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z \
46 +.083333333333333)/x
47
49 """Compute an incomplete gamma ratio.
50
51 Comments from Z. Yang:
52 Returns the incomplete gamma ratio I(x,alpha) where x is the upper
53 limit of the integration and alpha is the shape parameter.
54 returns (-1) if in error
55 ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
56 (1) series expansion if alpha>x or x<=1
57 (2) continued fraction otherwise
58 RATNEST FORTRAN by
59 Bhattacharjee GP (1970) The incomplete gamma integral. Applied Statistics,
60 19: 285-287 (AS32)
61 """
62 p = alpha
63 g = _ln_gamma_function(alpha)
64 accurate = 1e-8
65 overflow = 1e30
66 gin = 0
67 rn = 0
68 a = 0
69 b = 0
70 an = 0
71 dif = 0
72 term = 0
73 if x == 0:
74 return 0
75 if x < 0 or p <= 0:
76 return -1
77 factor = exp(p*log(x)-x-g)
78 if x > 1 and x >= p:
79 a = 1 - p
80 b = a + x + 1
81 term = 0
82 pn = [1, x, x+1, x*b, None, None]
83 gin = pn[2] / pn[3]
84 else:
85 gin=1
86 term=1
87 rn=p
88 while term > accurate:
89 rn += 1
90 term *= x / rn
91 gin += term
92 gin *= factor / p
93 return gin
94 while True:
95 a += 1
96 b += 2
97 term += 1
98 an = a * term
99 for i in range(2):
100 pn[i + 4] = b * pn[i + 2] - an * pn[i]
101 if pn[5] != 0:
102 rn = pn[4] / pn[5]
103 dif = abs(gin - rn)
104 if dif > accurate:
105 gin=rn
106 elif dif <= accurate*rn:
107 break
108 for i in range(4):
109 pn[i] = pn[i+2]
110 if abs(pn[4]) < overflow:
111 continue
112 for i in range(4):
113 pn[i] /= overflow
114 gin = 1 - factor * gin
115 return gin
116