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
|
module Distribution
module ChiSquare
module Ruby_
class << self
include Math
def pdf(x, n)
if n == 1
1.0 / Math.sqrt(2 * Math::PI * x) * Math::E**(-x / 2.0)
elsif n == 2
0.5 * Math::E**(-x / 2.0)
else
n = n.to_f
n2 = n / 2
x = x.to_f
1.0 / 2**n2 / gamma(n2) * x**(n2 - 1.0) * Math.exp(-x / 2.0)
end
end
# CDF Inverse over [x, \infty)
# Pr([x, \infty)) = y -> x
def pchi2(n, y)
if n == 1
w = Distribution::Normal.p_value(1 - y / 2) # = p1.0-Distribution::Normal.cdf(y/2)
w * w
elsif n == 2
# v = (1.0 / y - 1.0) / 33.0
# newton_a(y, v) {|x| [q_chi2(n, x), -chi2dens(n, x)] }
-2.0 * Math.log(y)
else
eps = 1.0e-5
v = 0.0
s = 10.0
loop do
v += s
break if s <= eps
if (qe = q_chi2(n, v) - y) == 0.0 then break end
if qe < 0.0
v -= s
s /= 10.0
end
end
v
end
end
def cdf(x, k)
1.0 - q_chi2(k, x)
end
# chi-square distribution ([1])
# Integral over [x, \infty)
def q_chi2(df, chi2)
chi2 = chi2.to_f
if (df & 1) != 0
chi = Math.sqrt(chi2)
return 2 * (1.0 - Distribution::Normal.cdf(chi)) if (df == 1)
s = t = chi * Math.exp(-0.5 * chi2) / SQ2PI
k = 3
while k < df
t *= chi2 / k; s += t
k += 2
end
2 * (1.0 - (Distribution::Normal.cdf(chi)) + s)
else
s = t = Math.exp(-0.5 * chi2)
k = 2
while k < df
t *= chi2 / k; s += t
k += 2
end
s
end
end
def quantile(pr, k)
pchi2(k, 1.0 - pr)
end
alias_method :p_value, :quantile
end
end
end
end
|
