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
|
# File 'lib/eps/statistics.rb', line 21
def self.incomplete_beta_function(x, alp, bet)
return if x < 0.0
return 1.0 if x > 1.0
tiny = 1.0E-50
if x > ((alp + 1.0)/(alp + bet + 2.0))
return 1.0 - incomplete_beta_function(1.0 - x, bet, alp)
end
lbet_ab = (Math.lgamma(alp)[0] + Math.lgamma(bet)[0] - Math.lgamma(alp + bet)[0]).freeze
front = (Math.exp(Math.log(x) * alp + Math.log(1.0 - x) * bet - lbet_ab) / alp.to_f).freeze
f, c, d = 1.0, 1.0, 0.0
returned_value = nil
(0..500).each do |number|
m = number/2
numerator = if number == 0
1.0
elsif number % 2 == 0
(m * (bet - m) * x)/((alp + 2.0 * m - 1.0)* (alp + 2.0 * m))
else
top = -((alp + m) * (alp + bet + m) * x)
down = ((alp + 2.0 * m) * (alp + 2.0 * m + 1.0))
top/down
end
d = 1.0 + numerator * d
d = tiny if d.abs < tiny
d = 1.0 / d
c = 1.0 + numerator / c
c = tiny if c.abs < tiny
cd = (c*d).freeze
f = f * cd
if (1.0 - cd).abs < 1.0E-10
returned_value = front * (f - 1.0)
break
end
end
returned_value
end
|