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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
module Math
def self.factorial(n)
return if n < 0
n = n.to_i # Only integers.
return 1 if n == 0 || n == 1
Math.gamma(n + 1) # Math.gamma(x) == (n - 1)! for integer values
end
def self.combination(n, r)
self.factorial(n)/(self.factorial(r) * self.factorial(n - r)).to_r # n!/(r! * [n - r]!)
end
def self.permutation(n, k)
self.factorial(n)/self.factorial(n - k).to_r
end
# Function adapted from the python implementation that exists in https://en.wikipedia.org/wiki/Simpson%27s_rule#Sample_implementation
# Finite integral in the interval [a, b] split up in n-intervals
def self.simpson_rule(a, b, n, &block)
unless n.even?
puts "The composite simpson's rule needs even intervals!"
return
end
h = (b - a)/n.to_r
resA = yield(a)
resB = yield(b)
sum = resA + resB
(1..n).step(2).each do |number|
res = yield(a + number * h)
sum += 4 * res
end
(1..(n-1)).step(2).each do |number|
res = yield(a + number * h)
sum += 2 * res
end
return sum * h / 3.0
end
def self.lower_incomplete_gamma_function(s, x)
base_iterator = x.round(1)
base_iterator += 1 if x < 1.0 && !x.zero?
# The greater the iterations, the better. That's why we are iterating 10_000 * x times
iterator = (10_000 * base_iterator).round
iterator = 100_000 if iterator.zero?
self.simpson_rule(0, x.to_r, iterator) do |t|
(t ** (s - 1)) * Math.exp(-t)
end
end
# Algorithm implementation translated from the ASA147 C++ version https://people.sc.fsu.edu/~jburkardt/cpp_src/asa147/asa147.html
# translated from FORTRAN by John Burkardt. Original algorithm written by Chi Leung Lau.
# It contains a modification on the error and underflow parameters to use maximum available float number
# and it performs the series using `Rational` objects to avoid memory exhaustion and reducing precision errors.
#
# This algorithm is licensed with MIT license.
#
# Reference:
#
# Chi Leung Lau,
# Algorithm AS 147:
# A Simple Series for the Incomplete Gamma Integral,
# Applied Statistics,
# Volume 29, Number 1, 1980, pages 113-114.
def self.normalised_lower_incomplete_gamma_function(s, x)
return 0.0 if s.negative? || x.zero? || x.negative?
# e = 1.0e-308
# uflo = 1.0e-47
e = Float::MIN
uflo = Float::MIN
lgamma, sign = Math.lgamma(s + 1.0)
arg = s * Math.log(x) - (sign * lgamma) - x
return 0.0 if arg < Math.log(uflo)
f = Math.exp(arg).to_r
return 0.0 if f.zero?
c = 1r
value = 1r
a = s.to_r
rational_x = x.to_r
loop do
a += 1r
c = c * (rational_x / a)
value += c
break if c <= (e * value)
end
(value * f).to_f
end
def self.beta_function(x, y)
return 1 if x == 1 && y == 1
(Math.gamma(x) * Math.gamma(y))/Math.gamma(x + y)
end
### This implementation is an adaptation of the incomplete beta function made in C by
### Lewis Van Winkle, which released the code under the zlib license.
### The whole math behind this code is described in the following post: https://codeplea.com/incomplete-beta-function-c
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 - self.incomplete_beta_function(1.0 - x, bet, alp)
end
# To avoid overflow problems, the implementation applies the logarithm properties
# to calculate in a faster and safer way the values.
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_r).freeze
# This is the non-log version of the left part of the formula (before the continuous fraction)
# down_left = alp * self.beta_function(alp, bet)
# upper_left = (x ** alp) * ((1.0 - x) ** bet)
# front = upper_left/down_left
f, c, d = 1.0, 1.0, 0.0
returned_value = nil
# Let's do more iterations than the proposed implementation (200 iters)
(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.to_r
c = 1.0 + numerator / c.to_r
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
end
|
