1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
"""Python implementation of chisqprob, to avoid SciPy dependency.
Adapted from SciPy: scipy/special/cephes/{chdtr,igam}.
"""
import math
# Cephes Math Library Release 2.0: April, 1987
# Copyright 1985, 1987 by Stephen L. Moshier
# Direct inquiries to 30 Frost Street, Cambridge, MA 02140
MACHEP = 0.0000001 # the machine roundoff error / tolerance
BIG = 4.503599627370496e15
BIGINV = 2.22044604925031308085e-16
def chisqprob(x, df):
"""Probability value (1-tail) for the Chi^2 probability distribution.
Broadcasting rules apply.
Parameters
----------
x : array_like or float > 0
df : array_like or float, probably int >= 1
Returns
-------
chisqprob : ndarray
The area from ``chisq`` to infinity under the Chi^2 probability
distribution with degrees of freedom ``df``.
"""
if x <= 0:
return 1.0
if x == 0:
return 0.0
if df <= 0:
raise ValueError("Domain error.")
if x < 1.0 or x < df:
return 1.0 - _igam(0.5 * df, 0.5 * x)
return _igamc(0.5 * df, 0.5 * x)
def _igamc(a, x):
"""Complemented incomplete Gamma integral (PRIVATE).
Parameters
----------
a: float
x: float
Returns
-------
float
Notes
-----
The function is defined by::
igamc(a,x) = 1 - igam(a,x)
inf.
-
1 | | -t a-1
= ----- | e t dt.
- | |
| (a) -
x
In this implementation both arguments must be positive.
The integral is evaluated by either a power series or
continued fraction expansion, depending on the relative
values of a and x.
"""
# Compute x**a * exp(-x) / Gamma(a)
ax = math.exp(a * math.log(x) - x - math.lgamma(a))
# Continued fraction
y = 1.0 - a
z = x + y + 1.0
c = 0.0
pkm2 = 1.0
qkm2 = x
pkm1 = x + 1.0
qkm1 = z * x
ans = pkm1 / qkm1
while True:
c += 1.0
y += 1.0
z += 2.0
yc = y * c
pk = pkm1 * z - pkm2 * yc
qk = qkm1 * z - qkm2 * yc
if qk != 0:
r = pk / qk
t = abs((ans - r) / r)
ans = r
else:
t = 1.0
pkm2 = pkm1
pkm1 = pk
qkm2 = qkm1
qkm1 = qk
if abs(pk) > BIG:
pkm2 *= BIGINV
pkm1 *= BIGINV
qkm2 *= BIGINV
qkm1 *= BIGINV
if t <= MACHEP:
return ans * ax
def _igam(a, x):
"""Left tail of incomplete Gamma function (PRIVATE).
Computes this formula::
inf. k
a -x - x
x e > ----------
- -
k=0 | (a+k+1)
"""
# Compute x**a * exp(-x) / Gamma(a)
ax = math.exp(a * math.log(x) - x - math.lgamma(a))
# Power series
r = a
c = 1.0
ans = 1.0
while True:
r += 1.0
c *= x / r
ans += c
if c / ans <= MACHEP:
return ans * ax / a
# --- Speed ---
# try:
# from scipy.stats import chisqprob
# except ImportError:
# pass
|
