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