Module: Distribution::MathExtension

Included in:
CMath, F::Ruby_, Math
Defined in:
lib/distribution/math_extension.rb,
lib/distribution/math_extension/erfc.rb,
lib/distribution/math_extension/gammastar.rb,
lib/distribution/math_extension/gsl_utilities.rb,
lib/distribution/math_extension/log_utilities.rb,
lib/distribution/math_extension/incomplete_beta.rb,
lib/distribution/math_extension/chebyshev_series.rb,
lib/distribution/math_extension/incomplete_gamma.rb,
lib/distribution/math_extension/exponential_integral.rb

Overview

Useful additions to Math

Defined Under Namespace

Modules: ApproxFactorial, Beta, Erfc, ExponentialIntegral, Gammastar, IncompleteBeta, IncompleteGamma, Log Classes: ChebyshevSeries, SwingFactorial

Constant Summary collapse

LNPI = 
1.14472988584940017414342735135
LN2 = 
0.69314718055994530941723212146
SQRT2 = 
1.41421356237309504880168872421
SQRTPI = 
1.77245385090551602729816748334
ROOT3_FLOAT_MIN = 
Float::MIN**(1 / 3.0)
ROOT3_FLOAT_EPSILON = 
Float::EPSILON**(1 / 3.0)
ROOT4_FLOAT_MIN = 
Float::MIN**(1 / 4.0)
ROOT4_FLOAT_EPSILON = 
Float::EPSILON**(1 / 4.0)
ROOT5_FLOAT_MIN = 
Float::MIN**(1 / 5.0)
ROOT5_FLOAT_EPSILON = 
Float::EPSILON**(1 / 5.0)
ROOT6_FLOAT_MIN = 
Float::MIN**(1 / 6.0)
ROOT6_FLOAT_EPSILON = 
Float::EPSILON**(1 / 6.0)
LOG_FLOAT_MIN = 
Math.log(Float::MIN)
EULER = 
0.57721566490153286060651209008

Instance Method Summary collapse

Instance Method Details

#beta(x, y) ⇒ Object

Beta function. Source:



177
178
179
 
# File 'lib/distribution/math_extension.rb', line 177

def beta(x, y)
  (gamma(x) * gamma(y)).quo(gamma(x + y))
end

#binomial_coefficient(n, k) ⇒ Object Also known as: combinations

Binomial coeffients, or: ( n ) ( k )

Gives the number of different k size subsets of a set size n

Uses:

(n) n^k' (n)..(n-k+1) ( ) = ---- = ------------ (k) k! k!



272
273
274
275
276
277
278
279
280
 
# File 'lib/distribution/math_extension.rb', line 272

def binomial_coefficient(n, k)
  return 1 if k == 0 || k == n
  k = [k, n - k].min
  permutations(n, k).quo(factorial(k))
  # The factorial way is
  # factorial(n).quo(factorial(k)*(factorial(n-k)))
  # The multiplicative way is
  # (1..k).inject(1) {|ac, i| (ac*(n-k+i).quo(i))}
end

#binomial_coefficient_gamma(n, k) ⇒ Object

Approximate binomial coefficient, using gamma function. The fastest method, until we fall on BigDecimal!



292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
 
# File 'lib/distribution/math_extension.rb', line 292

def binomial_coefficient_gamma(n, k)
  return 1 if k == 0 || k == n
  k = [k, n - k].min
  # First, we try direct gamma calculation for max precission

  val = gamma(n + 1).quo(gamma(k + 1) * gamma(n - k + 1))
  # Ups. Outside float point range. We try with logs
  if val.nan?
    # puts "nan"
    lg = loggamma(n + 1) - (loggamma(k + 1) + loggamma(n - k + 1))
    val = Math.exp(lg)
    # Crash again! We require BigDecimals
    if val.infinite?
      # puts "infinite"
      val = BigMath.exp(BigDecimal(lg.to_s), 16)
    end
  end
  val
end

#binomial_coefficient_multiplicative(n, k) ⇒ Object

Binomial coefficient using multiplicative algorithm On benchmarks, is faster that raising factorial method when k is little. Use only when you're sure of that.



284
285
286
287
288
 
# File 'lib/distribution/math_extension.rb', line 284

def binomial_coefficient_multiplicative(n, k)
  return 1 if k == 0 || k == n
  k = [k, n - k].min
  (1..k).inject(1) { |ac, i| (ac * (n - k + i).quo(i)) }
end

#erfc_e(x, with_error = false) ⇒ Object

Not the same as erfc. This is the GSL version, which may have slightly different results.



246
247
248
 
# File 'lib/distribution/math_extension.rb', line 246

def erfc_e(x, with_error = false)
  Erfc.evaluate(x, with_error)
end

#exact_regularized_beta(x, a, b) ⇒ Object

I_x(a,b): Regularized incomplete beta function TODO: Find a faster version. Source:



203
204
205
206
207
208
209
210
211
212
 
# File 'lib/distribution/math_extension.rb', line 203

def exact_regularized_beta(x, a, b)
  return 1 if x == 1

  m = a.to_i
  n = (b + a - 1).to_i

  (m..n).inject(0) do|sum, j|
    sum + (binomial_coefficient(n, j) * x**j * (1 - x)**(n - j))
  end
end

#exp_err(x, dx) ⇒ Object

e^x taking into account the error thus far (I think) gsl_sf_exp_err_e



24
25
26
27
28
29
 
# File 'lib/distribution/math_extension/gsl_utilities.rb', line 24

def exp_err(x, dx)
  adx = dx.abs
  fail('Overflow Error in exp_err: x + adx > LOG_FLOAT_MAX') if x + adx > LOG_FLOAT_MAX
  fail('Underflow Error in exp_err: x - adx < LOG_FLOAT_MIN') if x - adx < LOG_FLOAT_MIN
  [Math.exp(x), Math.exp(x) * [Float::EPSILON, Math.exp(adx) - 1.0 / Math.exp(adx)] + 2.0 * Float::EPSILON * Math.exp(x).abs]
end

#factorial(n) ⇒ Object

Exact factorial. Use lookup on a Hash table on n<20 and Prime Swing algorithm for higher values.



164
165
166
 
# File 'lib/distribution/math_extension.rb', line 164

def factorial(n)
  SwingFactorial.new(n).result
end

#fast_factorial(n) ⇒ Object

Approximate factorial, up to 16 digits Based of Luschy algorithm



170
171
172
 
# File 'lib/distribution/math_extension.rb', line 170

def fast_factorial(n)
  ApproxFactorial.stieltjes_factorial(n)
end

#gammq(a, x, with_error = false) ⇒ Object 



237
238
239
 
# File 'lib/distribution/math_extension.rb', line 237

def gammq(a, x, with_error = false)
  IncompleteGamma.q(a, x, with_error)
end

#incomplete_beta(x, a, b) ⇒ Object

Incomplete beta function: B(x;a,b) +a+ and +b+ are parameters and +x+ is integration upper limit.



217
218
219
220
 
# File 'lib/distribution/math_extension.rb', line 217

def incomplete_beta(x, a, b)
  IncompleteBeta.evaluate(a, b, x) * beta(a, b)
  # Math::IncompleteBeta.axpy(1.0, 0.0, a,b,x)
end

#incomplete_gamma(a, x = 0, with_error = false) ⇒ Object Also known as: gammp 



232
233
234
 
# File 'lib/distribution/math_extension.rb', line 232

def incomplete_gamma(a, x = 0, with_error = false)
  IncompleteGamma.p(a, x, with_error)
end

#lbeta(x, y) ⇒ Object

Log beta function conforming to style of lgamma (returns sign in second array index)



187
188
189
 
# File 'lib/distribution/math_extension.rb', line 187

def lbeta(x, y)
  Beta.log_beta(x, y)
end

#logbeta(x, y) ⇒ Object

Get pure-Ruby logbeta



182
183
184
 
# File 'lib/distribution/math_extension.rb', line 182

def logbeta(x, y)
  Beta.log_beta(x, y).first
end

#loggamma(x) ⇒ Object

Ln of gamma



228
229
230
 
# File 'lib/distribution/math_extension.rb', line 228

def loggamma(x)
  Math.lgamma(x).first
end

#permutations(n, k) ⇒ Object

Sequences without repetition. n^k' Also called 'failing factorial'



252
253
254
255
256
257
258
 
# File 'lib/distribution/math_extension.rb', line 252

def permutations(n, k)
  return 1 if k == 0
  return n if k == 1
  return factorial(n) if k == n
  (((n - k + 1)..n).inject(1) { |ac, v| ac * v })
  # factorial(x).quo(factorial(x-n))
end

#regularized_beta(x, a, b) ⇒ Object

I_x(a,b): Regularized incomplete beta function Fast version. For a exact calculation, based on factorial use exact_regularized_beta_function



194
195
196
197
 
# File 'lib/distribution/math_extension.rb', line 194

def regularized_beta(x, a, b)
  return 1 if x == 1
  IncompleteBeta.evaluate(a, b, x)
end

#rising_factorial(x, n) ⇒ Object

Rising factorial



223
224
225
 
# File 'lib/distribution/math_extension.rb', line 223

def rising_factorial(x, n)
  factorial(x + n - 1).quo(factorial(x - 1))
end

#unnormalized_incomplete_gamma(a, x, with_error = false) ⇒ Object 



241
242
243
 
# File 'lib/distribution/math_extension.rb', line 241

def unnormalized_incomplete_gamma(a, x, with_error = false)
  IncompleteGamma.unnormalized(a, x, with_error)
end