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
|
import math
import numpy as np
@np.vectorize
def norm_cdf(x: float, loc: float = 0.0, scale: float = 1.0) -> float:
x = (x - loc) / scale
x = x / 2**0.5
z = abs(x)
if z < 1 / 2**0.5:
y = 0.5 + 0.5 * math.erf(x)
else:
y = 0.5 * math.erfc(z)
if x > 0:
y = 1.0 - y
return y
@np.vectorize
def chi2_ppf(q: float) -> float:
"""
only deal with the special case df=1, loc=0, scale=1
solve chi2.cdf(x; df=1) = erf(sqrt(x/2)) = q with bisection method
"""
if q == 0:
return 0.0
if q == 1:
return math.inf
a, b = 0.0, 100.0
if q < 0.9:
for _ in range(100):
m = (a + b) / 2
if math.erf(math.sqrt(m / 2)) < q:
a = m
else:
b = m
else:
for _ in range(100):
m = (a + b) / 2
if math.erfc(math.sqrt(m / 2)) > 1.0 - q:
a = m
else:
b = m
return m
|
