Module: Distribution::MathExtension
- 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
-
#beta(x, y) ⇒ Object
Beta function.
-
#binomial_coefficient(n, k) ⇒ Object
(also: #combinations)
Binomial coeffients, or: ( n ) ( k ).
-
#binomial_coefficient_gamma(n, k) ⇒ Object
Approximate binomial coefficient, using gamma function.
-
#binomial_coefficient_multiplicative(n, k) ⇒ Object
Binomial coefficient using multiplicative algorithm On benchmarks, is faster that raising factorial method when k is little.
-
#erfc_e(x, with_error = false) ⇒ Object
Not the same as erfc.
-
#exact_regularized_beta(x, a, b) ⇒ Object
I_x(a,b): Regularized incomplete beta function TODO: Find a faster version.
-
#exp_err(x, dx) ⇒ Object
e^x taking into account the error thus far (I think) gsl_sf_exp_err_e.
-
#factorial(n) ⇒ Object
Exact factorial.
-
#fast_factorial(n) ⇒ Object
Approximate factorial, up to 16 digits Based of Luschy algorithm.
- #gammq(a, x, with_error = false) ⇒ Object
-
#incomplete_beta(x, a, b) ⇒ Object
Incomplete beta function: B(x;a,b)
a
andb
are parameters andx
is integration upper limit. - #incomplete_gamma(a, x = 0, with_error = false) ⇒ Object (also: #gammp)
-
#lbeta(x, y) ⇒ Object
Log beta function conforming to style of lgamma (returns sign in second array index).
-
#logbeta(x, y) ⇒ Object
Get pure-Ruby logbeta.
-
#loggamma(x) ⇒ Object
Ln of gamma.
-
#permutations(n, k) ⇒ Object
Sequences without repetition.
-
#regularized_beta(x, a, b) ⇒ Object
I_x(a,b): Regularized incomplete beta function Fast version.
-
#rising_factorial(x, n) ⇒ Object
Rising factorial.
- #unnormalized_incomplete_gamma(a, x, with_error = false) ⇒ Object
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 |