Module: Math

Defined in:
lib/math.rb

Class Method Summary collapse

Class Method Details

.beta_function(x, y) ⇒ Object 



107
108
109
110
111
 
# File 'lib/math.rb', line 107

def self.beta_function(x, y)
  return 1 if x == 1 && y == 1

  (Math.gamma(x) * Math.gamma(y))/Math.gamma(x + y)
end

.combination(n, r) ⇒ Object 



11
12
13
 
# File 'lib/math.rb', line 11

def self.combination(n, r)
  self.factorial(n)/(self.factorial(r) * self.factorial(n - r)).to_r # n!/(r! * [n - r]!)
end

.factorial(n) ⇒ Object 



2
3
4
5
6
7
8
9
 
# File 'lib/math.rb', line 2

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

.incomplete_beta_function(x, alp, bet) ⇒ Object

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: codeplea.com/incomplete-beta-function-c



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
 
# File 'lib/math.rb', line 116

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

.lower_incomplete_gamma_function(s, x) ⇒ Object 



47
48
49
50
51
52
53
54
55
56
57
58
 
# File 'lib/math.rb', line 47

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

.normalised_lower_incomplete_gamma_function(s, x) ⇒ Object

Algorithm implementation translated from the ASA147 C++ version 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.



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
 
# File 'lib/math.rb', line 74

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

.permutation(n, k) ⇒ Object 



15
16
17
 
# File 'lib/math.rb', line 15

def self.permutation(n, k)
  self.factorial(n)/self.factorial(n - k).to_r
end

.simpson_rule(a, b, n, &block) ⇒ Object

Function adapted from the python implementation that exists in en.wikipedia.org/wiki/Simpson%27s_rule#Sample_implementation Finite integral in the interval [a, b] split up in n-intervals



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
 
# File 'lib/math.rb', line 21

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