Module: Distribution::MathExtension

Included in:
CMath, Math
Defined in:
lib/distribution/math_extension.rb

Overview

Useful additions to Math

Defined Under Namespace

Modules: ApproxFactorial Classes: SwingFactorial

Instance Method Summary collapse

Instance Method Details

#beta(x, y) ⇒ Object

Beta function. Source:



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

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!



264
265
266
267
268
269
270
271
272
 
# File 'lib/distribution/math_extension.rb', line 264

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!



284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
 
# File 'lib/distribution/math_extension.rb', line 284

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.



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

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

#factorial(n) ⇒ Object

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



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

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

#fast_factorial(n) ⇒ Object

Approximate factorial, up to 16 digits Based of Luschy algorithm



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

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

#incomplete_beta(x, a, b) ⇒ Object

B_x(a,b) : Incomplete beta function Should be replaced by lib.stat.cmu.edu/apstat/63



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

def incomplete_beta(x,a,b)
  raise "Doesn't work"
end

#logbeta(x, y) ⇒ Object 



206
207
208
 
# File 'lib/distribution/math_extension.rb', line 206

def logbeta(x,y)
  (loggamma(x)+loggamma(y))-loggamma(x+y)
end

#loggamma(x) ⇒ Object

Ln of gamma



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

def loggamma(x)
  lg=Math.lgamma(x)
  lg[0]*lg[1]
end

#permutations(n, k) ⇒ Object

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



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

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_function(x, a, b) ⇒ Object

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



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

def regularized_beta_function(x,a,b)
  return 1 if x==1
  #incomplete_beta(x,a,b).quo(beta(a,b))
  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

#rising_factorial(x, n) ⇒ Object

Rising factorial



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

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