Module: Distribution::MathExtension

Included in:
CMath, 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:



220
221
222
 
# File 'lib/distribution/math_extension.rb', line 220

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!



313
314
315
316
317
318
319
320
321
 
# File 'lib/distribution/math_extension.rb', line 313

def binomial_coefficient(n,k)
  return 1 if (k==0 or 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!



333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
 
# File 'lib/distribution/math_extension.rb', line 333

def binomial_coefficient_gamma(n,k)
  return 1 if (k==0 or 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.



325
326
327
328
329
 
# File 'lib/distribution/math_extension.rb', line 325

def binomial_coefficient_multiplicative(n,k)
  return 1 if (k==0 or 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.



287
288
289
 
# File 'lib/distribution/math_extension.rb', line 287

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:



245
246
247
248
249
250
251
252
253
 
# File 'lib/distribution/math_extension.rb', line 245

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) {|sum,j|
   sum+(binomial_coefficient(n,j)* x**j * (1-x)**(n-j))
 }

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.



208
209
210
 
# File 'lib/distribution/math_extension.rb', line 208

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

#fast_factorial(n) ⇒ Object

Approximate factorial, up to 16 digits Based of Luschy algorithm



213
214
215
 
# File 'lib/distribution/math_extension.rb', line 213

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

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



278
279
280
 
# File 'lib/distribution/math_extension.rb', line 278

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.



258
259
260
261
 
# File 'lib/distribution/math_extension.rb', line 258

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 



273
274
275
 
# File 'lib/distribution/math_extension.rb', line 273

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)



230
231
232
 
# File 'lib/distribution/math_extension.rb', line 230

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

#logbeta(x, y) ⇒ Object

Get pure-Ruby logbeta



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

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

#loggamma(x) ⇒ Object

Ln of gamma



269
270
271
 
# File 'lib/distribution/math_extension.rb', line 269

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

#permutations(n, k) ⇒ Object

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



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

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



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

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

#rising_factorial(x, n) ⇒ Object

Rising factorial



264
265
266
 
# File 'lib/distribution/math_extension.rb', line 264

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

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



282
283
284
 
# File 'lib/distribution/math_extension.rb', line 282

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