|
1
|
|
|
# Copyright (c) 2014, Salesforce.com, Inc. All rights reserved. |
|
2
|
|
|
# Copyright (c) 2015-2016, Gamelan Labs, Inc. |
|
3
|
|
|
# Copyright (c) 2016, Google, Inc. |
|
4
|
|
|
# |
|
5
|
|
|
# Redistribution and use in source and binary forms, with or without |
|
6
|
|
|
# modification, are permitted provided that the following conditions |
|
7
|
|
|
# are met: |
|
8
|
|
|
# |
|
9
|
|
|
# - Redistributions of source code must retain the above copyright |
|
10
|
|
|
# notice, this list of conditions and the following disclaimer. |
|
11
|
|
|
# - Redistributions in binary form must reproduce the above copyright |
|
12
|
|
|
# notice, this list of conditions and the following disclaimer in the |
|
13
|
|
|
# documentation and/or other materials provided with the distribution. |
|
14
|
|
|
# - Neither the name of Salesforce.com nor the names of its contributors |
|
15
|
|
|
# may be used to endorse or promote products derived from this |
|
16
|
|
|
# software without specific prior written permission. |
|
17
|
|
|
# |
|
18
|
|
|
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
|
19
|
|
|
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
|
20
|
|
|
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
|
21
|
|
|
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE |
|
22
|
|
|
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, |
|
23
|
|
|
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, |
|
24
|
|
|
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS |
|
25
|
|
|
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND |
|
26
|
|
|
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR |
|
27
|
|
|
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE |
|
28
|
|
|
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
|
29
|
|
|
|
|
30
|
|
|
""" |
|
31
|
|
|
This module computes the chi^2 survival function and the required |
|
32
|
|
|
functions. |
|
33
|
|
|
""" |
|
34
|
|
|
|
|
35
|
|
|
from __future__ import division |
|
36
|
|
|
import math |
|
37
|
|
|
import sys |
|
38
|
|
|
|
|
39
|
|
|
|
|
40
|
|
|
def log(x): |
|
41
|
|
|
if x > sys.float_info.min: |
|
42
|
|
|
value = math.log(x) |
|
43
|
|
|
else: |
|
44
|
|
|
value = -float('inf') |
|
45
|
|
|
return value |
|
46
|
|
|
|
|
47
|
|
|
|
|
48
|
|
|
def incomplete_gamma(x, s): |
|
49
|
|
|
r""" |
|
50
|
|
|
This function computes the incomplete lower gamma function |
|
51
|
|
|
using the series expansion: |
|
52
|
|
|
|
|
53
|
|
|
.. math:: |
|
54
|
|
|
|
|
55
|
|
|
\gamma(x, s) = x^s \Gamma(s)e^{-x}\sum^\infty_{k=0} |
|
56
|
|
|
\frac{x^k}{\Gamma(s + k + 1)} |
|
57
|
|
|
|
|
58
|
|
|
This series will converge strongly because the Gamma |
|
59
|
|
|
function grows factorially. |
|
60
|
|
|
|
|
61
|
|
|
Because the Gamma function does grow so quickly, we can |
|
62
|
|
|
run into numerical stability issues. To solve this we carry |
|
63
|
|
|
out as much math as possible in the log domain to reduce |
|
64
|
|
|
numerical error. This function matches the results from |
|
65
|
|
|
scipy to numerical precision. |
|
66
|
|
|
""" |
|
67
|
|
|
if x < 0: |
|
68
|
|
|
return 1 |
|
69
|
|
|
if x > 1e3: |
|
70
|
|
|
return math.gamma(s) |
|
71
|
|
|
log_gamma_s = math.lgamma(s) |
|
72
|
|
|
log_x = log(x) |
|
73
|
|
|
value = 0 |
|
74
|
|
|
for k in range(100): |
|
75
|
|
|
log_num = (k + s)*log_x + (-x) + log_gamma_s |
|
76
|
|
|
log_denom = math.lgamma(k + s + 1) |
|
77
|
|
|
value += math.exp(log_num - log_denom) |
|
78
|
|
|
return value |
|
79
|
|
|
|
|
80
|
|
|
|
|
81
|
|
|
def chi2sf(x, s): |
|
82
|
|
|
r""" |
|
83
|
|
|
This function returns the survival function of the chi^2 |
|
84
|
|
|
distribution. The survival function is given as: |
|
85
|
|
|
|
|
86
|
|
|
.. math:: |
|
87
|
|
|
1 - CDF |
|
88
|
|
|
|
|
89
|
|
|
where rhe chi^2 CDF is given as |
|
90
|
|
|
|
|
91
|
|
|
.. math:: |
|
92
|
|
|
F(x; s) = \frac{ \gamma( x/2, s/2 ) }{ \Gamma(s/2) }, |
|
93
|
|
|
|
|
94
|
|
|
with :math:`\gamma` is the incomplete gamma function defined above. |
|
95
|
|
|
Therefore, the survival probability is givne by: |
|
96
|
|
|
|
|
97
|
|
|
.. math:: |
|
98
|
|
|
1 - \frac{ \gamma( x/2, s/2 ) }{ \Gamma(s/2) }. |
|
99
|
|
|
|
|
100
|
|
|
This function matches the results from |
|
101
|
|
|
scipy to numerical precision. |
|
102
|
|
|
""" |
|
103
|
|
|
top = incomplete_gamma(x/2, s/2) |
|
104
|
|
|
bottom = math.gamma(s/2) |
|
105
|
|
|
value = top/bottom |
|
106
|
|
|
return 1 - value |
|
107
|
|
|
|
posterior /
goftests
