|
1
|
|
|
import math |
|
2
|
|
|
import sys |
|
3
|
|
|
from numpy import inf |
|
4
|
|
|
|
|
5
|
|
|
|
|
6
|
|
|
''' |
|
7
|
|
|
This was implemented so that we could use the goftests module |
|
8
|
|
|
with pypy. Given that the only scipy dependence inside goftests |
|
9
|
|
|
was from the chi^2 survival function (1 - cdf), and then gamma |
|
10
|
|
|
function. |
|
11
|
|
|
|
|
12
|
|
|
This module computes the chi^2 survival function and the requied |
|
13
|
|
|
functions. |
|
14
|
|
|
|
|
15
|
|
|
Copyright (c) 2016, Gamelan Labs, Inc. |
|
16
|
|
|
''' |
|
17
|
|
|
|
|
18
|
|
|
|
|
19
|
|
|
def log(x): |
|
20
|
|
|
if x > sys.float_info.min: |
|
21
|
|
|
value = math.log(x) |
|
22
|
|
|
else: |
|
23
|
|
|
value = -inf |
|
24
|
|
|
return value |
|
25
|
|
|
|
|
26
|
|
|
|
|
27
|
|
|
def incomplete_gamma(x, s): |
|
28
|
|
|
''' |
|
29
|
|
|
This function computes the incomplete lower gamma function |
|
30
|
|
|
using the series expansion: |
|
31
|
|
|
\gamma(x, s) = x^s \Gamma(s)e^{-x}\sum^\infty_{k=0} |
|
32
|
|
|
\frac{x^k}{\Gamma(s + k + 1)} |
|
33
|
|
|
This series will converge strongly because the Gamma |
|
34
|
|
|
function grows factorially. |
|
35
|
|
|
|
|
36
|
|
|
Because the Gamma function does grow so quickly, we can |
|
37
|
|
|
run into numerical stability issues. I have chosen to carry |
|
38
|
|
|
out as much math as possible in the log domain to reduce |
|
39
|
|
|
numerical error. This function matches the results from |
|
40
|
|
|
scipy to numerical precision. |
|
41
|
|
|
|
|
42
|
|
|
''' |
|
43
|
|
|
value = 0.0 |
|
44
|
|
|
if x < 0.0: |
|
45
|
|
|
return 1.0 |
|
46
|
|
|
if x > 1e3: |
|
47
|
|
|
return math.gamma(s) |
|
48
|
|
|
log_gamma_s = log(math.gamma(s)) |
|
49
|
|
|
log_x = log(x) |
|
50
|
|
|
for k in range(0, 100): |
|
51
|
|
|
log_num = (k + s)*log_x + (-x) + log_gamma_s |
|
52
|
|
|
log_denom = log(math.gamma(k + s + 1)) |
|
53
|
|
|
value += math.exp(log_num - log_denom) |
|
54
|
|
|
return value |
|
55
|
|
|
|
|
56
|
|
|
|
|
57
|
|
|
def chi2sf(x, s): |
|
58
|
|
|
''' |
|
59
|
|
|
This function returns the survival function of the chi^2 |
|
60
|
|
|
distribution. The survival function is given as: |
|
61
|
|
|
1 - CDF |
|
62
|
|
|
where rhe chi^2 CDF is given as |
|
63
|
|
|
F(x; s) = \frac{ \gamma( x/2, s/2 ) }{ \Gamma(s/2) }, |
|
64
|
|
|
with $\gamma$ is the incomplete gamma function defined above. |
|
65
|
|
|
Therefore, the survival probability is givne by: |
|
66
|
|
|
1 - \frac{ \gamma( x/2, s/2 ) }{ \Gamma(s/2) }. |
|
67
|
|
|
This function matches the results from |
|
68
|
|
|
scipy to numerical precision. |
|
69
|
|
|
''' |
|
70
|
|
|
top = incomplete_gamma(x/2.0, s/2.0) |
|
71
|
|
|
bottom = math.gamma(s/2.0) |
|
72
|
|
|
value = top/bottom |
|
73
|
|
|
return 1.0 - value |
|
74
|
|
|
|
posterior /
goftests
