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:



185
186
187
 
# File 'lib/distribution/math_extension.rb', line 185

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!



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

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!



301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
 
# File 'lib/distribution/math_extension.rb', line 301

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.



293
294
295
296
297
 
# File 'lib/distribution/math_extension.rb', line 293

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.



255
256
257
 
# File 'lib/distribution/math_extension.rb', line 255

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:



211
212
213
214
215
216
217
218
219
220
 
# File 'lib/distribution/math_extension.rb', line 211

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
  raise("Overflow Error in exp_err: x + adx > LOG_FLOAT_MAX") if x + adx > LOG_FLOAT_MAX
  raise("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.



172
173
174
 
# File 'lib/distribution/math_extension.rb', line 172

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

#fast_factorial(n) ⇒ Object

Approximate factorial, up to 16 digits Based of Luschy algorithm



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

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

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



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

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.



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

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 



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

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)



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

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

#logbeta(x, y) ⇒ Object

Get pure-Ruby logbeta



190
191
192
 
# File 'lib/distribution/math_extension.rb', line 190

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

#loggamma(x) ⇒ Object

Ln of gamma



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

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

#permutations(n, k) ⇒ Object

Sequences without repetition. n^k’ Also called ‘failing factorial’



261
262
263
264
265
266
267
 
# File 'lib/distribution/math_extension.rb', line 261

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



202
203
204
205
 
# File 'lib/distribution/math_extension.rb', line 202

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

#rising_factorial(x, n) ⇒ Object

Rising factorial



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

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

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



250
251
252
 
# File 'lib/distribution/math_extension.rb', line 250

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