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
|
module Distribution
module T
module Ruby_
class << self
def pdf(t,v)
((Math.gamma((v+1) / 2.0)) / (Math.sqrt(v*Math::PI)*Math.gamma(v/2.0))) * ((1+(t**2 / v.to_f))**(-(v+1) / 2.0))
end
# Returns the integral of t-distribution with n degrees of freedom over (-Infty, x].
def cdf(t, n)
p_t(n, t)
end
# t-distribution ([1])
# (-\infty, x]
def p_t(df, t)
if df.to_i!=df
x=(t+Math.sqrt(t**2+df)) / (2*Math.sqrt(t**2+df))
return Math.regularized_beta(x,df/2.0,df/2.0)
end
df=df.to_i
c2 = df.to_f / (df + t * t);
s = Math.sqrt(1.0 - c2)
s = -s if t < 0.0
p = 0.0;
i = df % 2 + 2
while i <= df
p += s
s *= (i - 1) * c2 / i
i += 2
end
if df.is_a? Float or df & 1 != 0
0.5+(p*Math.sqrt(c2)+Math.atan(t/Math.sqrt(df))) / Math::PI
else
(1.0 + p) / 2.0
end
end
# inverse of t-distribution ([2])
# (-\infty, -q/2] + [q/2, \infty)
def ptsub(q, n)
q = q.to_f
if(n == 1 && 0.001 < q && q < 0.01)
eps = 1.0e-4
elsif (n == 2 && q < 0.0001)
eps = 1.0e-4
elsif (n == 1 && q < 0.001)
eps = 1.0e-2
else
eps = 1.0e-5
end
s = 10000.0
w = 0.0
loop do
w += s
if(s <= eps) then return w end
if((qe = 2.0 - p_t(n, w)*2.0 - q) == 0.0) then return w end
if(qe < 0.0)
w -= s
s /= 10.0 #/
end
end
end
def pt(q, n)
q = q.to_f
if(q < 1.0e-5 || q > 1.0 || n < 1)
$stderr.printf("Error : Illigal parameter in pt()!\n")
return 0.0
end
if(n <= 5) then return ptsub(q, n) end
if(q <= 5.0e-3 && n <= 13) then return ptsub(q, n) end
f1 = 4.0 * (f = n.to_f)
f5 = (f4 = (f3 = (f2 = f * f) * f) * f) * f
f2 *= 96.0
f3 *= 384.0
f4 *= 92160.0
f5 *= 368640.0
u = Normal.p_value(1.0 - q / 2.0)
w0 = (u2 = u * u) * u
w1 = w0 * u2
w2 = w1 * u2
w3 = w2 * u2
w4 = w3 * u2
w = (w0 + u) / f1
w += (5.0 * w1 + 16.0 * w0 + 3.0 * u) / f2
w += (3.0 * w2 + 19.0 * w1 + 17.0 * w0 - 15.0 * u) / f3
w += (79.0 * w3 + 776.0 * w2 + 1482.0 * w1 - 1920.0 * w0 - 9450.0 * u) / f4
w += (27.0 * w4 + 339.0 * w3 + 930.0 * w2 - 1782.0 * w1 - 765.0 * w0 + 17955.0 * u) / f5
u + w
end
# Returns the P-value of tdist().
def p_value(y,n)
if y > 0.5
pt(2.0 - y*2.0, n)
else
- pt(y*2.0, n)
end
end
end
end
end
end
|
