# module Math

### Constants

EPSILON
FACTORIALS

  First 16 factorials.

INVERSE_LN_2

### Public Class Methods

abs(x) click to toggle source

Absolute value of `x`.

    # File lib/facets/math/abs.rb, line 4
    def self.abs(x)
      x.abs
    end

absolute_mean_difference(array)
Alias for: amd
acot(x) click to toggle source

Arcus cotangens of `x`

    # File lib/facets/math/acot.rb, line 4
    def self.acot(x)
      (PI * 0.5) - atan(x)
    end

acoth(x) click to toggle source

Area cotangens hyperbolicus of `x`

    # File lib/facets/math/acoth.rb, line 4
    def self.acoth(x)
      0.5 * log((x + 1.0) / (x - 1.0))
    end

acsc(x) click to toggle source

Arcus cosecans of `x`

    # File lib/facets/math/acsc.rb, line 4
    def self.acsc(x)
      asin(1.0 / x)
    end

amd(array) click to toggle source

The average absolute difference of two independent values drawn from the sample. Equal to the RMD * mean.

    # File lib/facets/math/amd.rb, line 8
    def self.amd(array)
            rmd(array) * mean(array)
    end

Also aliased as: absolute_mean_difference
approx_equal(a, b, epsilon=EPSILON) click to toggle source

Approximately equal.

TODO: Use core extension Numeric#approx? instead (?)

    # File lib/facets/math/approx_equal.rb, line 9
    def self.approx_equal(a, b, epsilon=EPSILON)
            c = a - b
            c *= -1.0 if c < 0
      c < epsilon
    end

asec(x) click to toggle source

Arcus secans of `x`

    # File lib/facets/math/asec.rb, line 4
    def self.asec(x)
      acos(1.0 / x)
    end

atkinson_index(array) click to toggle source

Closely related to the Theil index and easily expressible in terms of it.

AI = 1-e^{theil_index}

en.wikipedia.org/wiki/Atkinson_index

    # File lib/facets/math/atkinson_index.rb, line 11
    def self.atkinson_index(array)
      t = theil_index(array)
      (t < 0) ? -1 : 1-Math::E**(-t)
    end

beta(x, y) click to toggle source

Beta function of `x` and `y` - `beta(x, y) = tgamma(x) * tgamma(y) / tgamma(x + y)`

    # File lib/facets/math/beta.rb, line 5
    def self.beta(x, y)
      exp(lgamma(x) + lgamma(y) - lgamma(x+y))
    end

cdf(array, normalised=1.0) click to toggle source

Returns the Cumulative Density Function of this sample (normalised to a fraction of 1.0).

    # File lib/facets/math/cdf.rb, line 5
    def self.cdf(array, normalised=1.0)
      s = sum(array).to_f
      array.sort.inject([0.0]) { |c,d| c << c[-1] + normalised*d.to_f/s }
    end

ceil(x) click to toggle source

Smallest integer not smaller than `x`.

    # File lib/facets/math/ceil.rb, line 4
    def self.ceil(x)
      x.ceil
    end

cot(x) click to toggle source

Cotangens of `x`

    # File lib/facets/math/cot.rb, line 4
    def self.cot(x)
      tan((PI * 0.5) - x)
    end

coth(x) click to toggle source

Cotangens hyperbolicus of `x`

    # File lib/facets/math/coth.rb, line 4
    def self.coth(x)
      1.0 / tanh(x)
    end

csc(x) click to toggle source

Cosecans of `x`

    # File lib/facets/math/csc.rb, line 4
    def self.csc(x)
      1.0 / sin(x)
    end

csch(x) click to toggle source

Cosecans hyperbolicus of `x`

    # File lib/facets/math/csch.rb, line 4
    def self.csch(x)
      1.0 / sinh(x)
    end

delta(i, j) click to toggle source

Kronecker symbol of `i` and `j`. Returns 1 if `i` and `j` are equal, 0 otherwise.

    # File lib/facets/math/delta.rb, line 5
    def self.delta(i, j)
      return Integer(i) == Integer(j) ? 1 : 0
    end

epsilon(i, j, k) click to toggle source

Levi-Civita symbol of `i`, `j`, and `k` - 1 if (`i`, `j`, `k`) is (1, 2, 3), (2, 3, 1), or (3, 1, 2), -1 if it is (1, 3, 2), (2, 1, 3), or (3, 2, 1), 0 as long as `i`, `j`, and `k` are all elements of {1, 2, 3}, otherwise returns `nil`.

    # File lib/facets/math/epsilon.rb, line 7
    def self.epsilon(i, j, k)
      i = Integer(i)
      return nil if i < 1 or i > 3
      j = Integer(j)
      return nil if j < 1 or j > 3
      k = Integer(k)
      return nil if k < 1 or k > 3
      case i * 16 + j * 4 + k
        when 27, 45, 54 then return  1
        when 30, 39, 57 then return -1
      end
      0
    end

exp10(x) click to toggle source

10 to the power `x`

    # File lib/facets/math/exp10.rb, line 4
    def self.exp10(x)
      10.0 ** x
    end

exp2(x) click to toggle source

2 to the power `x`

    # File lib/facets/math/exp2.rb, line 4
    def self.exp2(x)
      2.0 ** x
    end

factorial(n) click to toggle source

1 * 2 * … * `n`, `nil` for negative numbers

    # File lib/facets/math/factorial.rb, line 24
    def self.factorial(n)
      n = Integer(n)
      if n < 0
        nil
      elsif FACTORIALS.length > n
        FACTORIALS[n]
      else
        h = FACTORIALS.last
        (FACTORIALS.length .. n).each { |i| FACTORIALS.push h *= i }
        h
      end
    end

floor(x) click to toggle source

Largest integer not larger than `x`.

    # File lib/facets/math/floor.rb, line 4
    def self.floor(x)
      x.floor
    end

gamma(x) click to toggle source

    # File lib/facets/math/gamma.rb, line 7
    def self.gamma(x)
      exp(lgamma(x))
    end

gcd(m, n) click to toggle source

Greatest common divisor of `m` and `n`, `nil` for non-positive numbers - gcd is computed by means of the Euclidian algorithm.

    # File lib/facets/math/gcd.rb, line 5
    def self.gcd(m, n)
      m = Integer(m)
      n = Integer(n)
      if m <= 0 || n <= 0
        return nil
      end
      loop {
        if m < n
          m, n = n, m
        end
        if (l = m % n) == 0
          break
        end
        m = l
      }
      n
    end

gini_coefficient(array) click to toggle source

Calculates the Gini Coefficient (a measure of inequality of a distribution based on the area between the Lorenz curve and the uniform curve).

en.wikipedia.org/wiki/Gini_coefficient

This is a slightly cleaner way of calculating the Gini Coefficient then the previous implementationj.

    GC = \frac{\sum_{i=1}^N (2i-N-1)x_i}{N^2-\bar{x}}

    # File lib/facets/math/gini_coefficient.rb, line 15
    def self.gini_coefficient(array)
      return -1 if size <= 0 or any? { |x| x < 0 }
      return 0 if size < 2 or all? { |x| approx_equal(x,0) }
      s = 0
      sort.each_with_index { |li,i| s += (2*i+1-size)*li }
      s.to_f/(size**2*mean).to_f
    end

kldivergence(array, q) click to toggle source

The Kullback-Leibler divergence from this array to that of `q`.

NB: You will possibly want to sort both P and Q before calling this depending on what you're actually trying to measure.

en.wikipedia.org/wiki/Kullback-Leibler_divergence

    # File lib/facets/math/kldivergence.rb, line 10
    def self.kldivergence(array, q)
      fail "Buggy."
      fail "Cannot compare differently sized arrays." unless size = q.size
      kld = 0
      each_with_index { |pi,i| kld += pi*Math::log(pi.to_f/q[i].to_f) }
      kld
    end

lcm(m, n) click to toggle source

Least common multiple of `m` and `n`, computed by multiplying `m` and `n` and dividing the product by the gcd of `m` and `n`, `nil` for non-positive numbers.

    # File lib/facets/math/lcm.rb, line 6
    def self.lcm(m, n)
      m = Integer(m)
      n = Integer(n)
      if m <= 0 || n <= 0
        return nil
      end
      m / gcd(m, n) * n
    end

lgamma(x) click to toggle source

Logarithmus naturalis of gamma function of `x`.

@author Josef Schugt

    # File lib/facets/math/lgamma.rb, line 8
    def self.lgamma(x)
      h  = x + 5.5
      h -= (x + 0.5) * log(h)
      sum  =  1.000_000_000_190_015
      sum += 76.180_091_729_471_46           / (x + 1.0)
      sum -= 86.505_320_329_416_77           / (x + 2.0)
      sum += 24.014_098_240_830_91           / (x + 3.0)
      sum -=  1.231_739_572_450_155          / (x + 4.0)
      sum +=  0.120_865_097_386_617_9e-2 / (x + 5.0)
      sum -=  0.539_523_938_495_3e-5     / (x + 6.0)
      -h + log(2.506_628_274_631_000_5 * sum / x)
    end

linsolve(a, b, c = 0.0) click to toggle source

Returns real solution(s) of `+a+x + b = c` or `nil` if no or an infinite number of solutions exist. If `c` is missing it is assumed to be 0.

@author Josef Schugt

    # File lib/facets/math/linsolve.rb, line 8
    def self.linsolve(a, b, c = 0.0)
      a == 0 ? nil : (c - b) / a
    end

log2(x) click to toggle source

Logarithmus dualis of `x`.

    # File lib/facets/math/log2.rb, line 8
    def self.log2(x)
      Math.log(x) * INVERSE_LN_2
    end

max(array, block) { |max| ... } click to toggle source

    # File lib/facets/math/min.rb, line 20
    def self.max(array, block)
      if block_given?
        if max = find{|i| i}
          max = yield(max)
          each{|i|
            j = yield(i)
            max = j if max < j
          }
          max
        end
      else
        array.max
      end
    end

mean(array, &blk) click to toggle source

Mean average.

    # File lib/facets/math/mean.rb, line 6
    def self.mean(array, &blk)
      s = array.size
      return 0.0 if s == 0
      sum(array, &blk) / s
    end

Also aliased as: mean_average
mean_average(array, &blk)
Alias for: mean
median(array) click to toggle source

Returns the numerical median for the an array of values; or nil if array is empty.

    # File lib/facets/math/median.rb, line 8
    def self.median(array)
      percentile(array, 50)
    end

min(array) { |min| ... } click to toggle source

    # File lib/facets/math/min.rb, line 4
    def self.min(array, &block)
      if block_given?
        if min = array.find{ |i| i }
          min = yield(min)
          array.each do |i|
            j = yield(i)
            min = j if min > j
          end
          min
        end
      else
        array.min
      end
    end

percentile(array, pcnt) click to toggle source

Returns the percentile value for percentile pcnt; nil if array is empty.

`pcnt` should be expressed as an integer, e.g. `percentile(90)` returns the 90th percentile of the array.

Algorithm from NIST

NOTE: This is not a common core extension and is not loaded automatically when using `require 'facets'`.

CREDIT: Ben Koski

@non-core

    require 'facets/array/precentile'

    # File lib/facets/math/percentile.rb, line 18
    def self.percentile(array, pcnt)
      sorted_array = array.sort

      return nil if array.length == 0

      rank  = (pcnt.to_f / 100) * (array.length + 1)
      whole = rank.truncate

      # if has fractional part
      if whole != rank
        s0 = sorted_array[whole - 1]
        s1 = sorted_array[whole]

        f = (rank - rank.truncate).abs

        return (f * (s1 - s0)) + s0
      else
        return sorted_array[whole - 1]
      end
    end

pow(x, y) click to toggle source

`x` to the power `y`

    # File lib/facets/math/pow.rb, line 4
    def self.pow(x, y)
      x ** y
    end

pstd(array, &block) click to toggle source

Standard deviation of a population.

    # File lib/facets/math/std.rb, line 17
    def self.pstd(array, &block)
      Math::sqrt(pvariance(array, &block))
    end

pvariance(array) click to toggle source

Variance of a population. Variance of 0 or 1 elements is 0.0.

    # File lib/facets/math/variance.rb, line 26
    def self.pvariance(array)
      return 0.0 if array.size < 2
      summed_sqdevs(array) / array.size
    end

relative_mean_difference(array)
Alias for: rmd
rmd(array) click to toggle source

Calculates the relative mean difference of this sample. Makes use of the fact that the Gini Coefficient is half the RMD.

    # File lib/facets/math/rmd.rb, line 7
    def self.rmd(array)
      return 0.0 if approx_equal(mean(array), 0.0)
      gini_coefficient(array) * 2
    end

Also aliased as: relative_mean_difference
root(x, y) click to toggle source

The `y` root of `x`.

    # File lib/facets/math/root.rb, line 4
    def self.root(x, y)
      x ** (1.0 / y)
    end

round(x) click to toggle source

Round number to an integer.

    # File lib/facets/math/round.rb, line 5
    def self.round(x)
      x.round
    end

sec(x) click to toggle source

Secans of `x`.

    # File lib/facets/math/sec.rb, line 4
    def self.sec(x)
      1.0 / cos(x)
    end

sech(x) click to toggle source

Secans hyperbolicus of `x`

    # File lib/facets/math/sech.rb, line 4
    def self.sech(x)
      1.0 / cosh(x)
    end

sign(x) click to toggle source

Sign of `x`. Returns -1 for negative x, +1 for positive x and zero for x = 0

    # File lib/facets/math/sign.rb, line 5
    def self.sign(x)
      (x > 0.0) ? 1.0 : ((x < 0.0) ? -1.0 : 0.0)
    end

sinc(x) click to toggle source

Sinc function of `x`.

    # File lib/facets/math/sinc.rb, line 4
    def self.sinc(x)
      (x == 0.0) ? 1.0 : sin(x) / x
    end

sqr(x) click to toggle source

Square of number.

    # File lib/facets/math/sqr.rb, line 4
    def self.sqr(x)
      x * x
    end

sqsolve(a, b, c, d = 0.0) click to toggle source

Returns array of real solution of `ax**2 + bx + c = d` or `nil` if no or an infinite number of solutions exist. If `d` is missing it is assumed to be 0.

In order to solve `ax**2 + bx + c = d` `sqsolve` identifies several cases:

- `a == 0:` The equation to be solved is the linear equation `bx + c = d`. sqsolve> delegates the computation to linsolve>. If it results in `nil`, `nil` is returned (not `[nil]`!). Otherwise a one-element array containing result of linsolve is returned.

- `a != 0:`

      The equation to be solved actually is a second order one.
      * <code>c == d</code>
        The equation to be solved is <code>ax**2 + bx = 0</code>. One solution of this equation obviously is
        <code>x = 0</code>, the second one solves <code>ax + b = 0</code>. The solution of the latter is
        delegated to +linsolve+. An array containing both results in ascending order is returned.
      * <code>c != d</code>
        The equation cannot be separated into <code>x</code> times some factor.
        * <code>b == 0</code>
          The equation to be solved is <code>ax**2 + c = d</code>. This can be written as the linear equation
          <code>ay + c = d</code> with <code>y = x ** 2</code>. The solution of the linear equation is delegated
          to +linsolve+. If the returned value for +y+ is +nil+, that becomes the overall return value.
          Otherwise an array containing the negative and positive squareroot of +y+ is returned
        * <code>b != 0 </code>
          The equation cannot be reduced to simpler cases. We now first have to compute what is called the
          discriminant <code>x = b**2 + 4a(d - c)</code> (that's what we need to compute the square root of).
          If the descriminant is negative no real solution exists and <code>nil</code> is returned. The ternary
          operator checking whether <code>b</code> is negative does ensure better numerical stability --only one
          of the two solutions is computed using the widely know formula for solving second order equations.
          The second one is computed from the fact that the product of both solutions is <code>(c - d) / a</code>.
          Take a look at a book on numerical mathematics if you don't understand why this should be done.

@author Josef Schugt

    # File lib/facets/math/sqsolve.rb, line 37
    def self.sqsolve(a, b, c, d = 0.0)
      if a == 0.0
        x = linsolve(b, c, d)
        return x.nil? ? nil: [ linsolve(b, c, d) ]
      else
        return [0.0, linsolve(a, b)].sort if c == d
        if b == 0.0
          x = linsolve(a, c, d)
          x < 0.0 ? nil : [-Math.sqrt(x), Math.sqrt(x)]
        else
          x = b * b + 4.0 * a * (d - c)
          return nil if x < 0.0
          x = b < 0 ? b - Math.sqrt(x) : b + Math.sqrt(x)
          [-0.5 * x / a, 2.0 * (d - c) / x].sort
        end
      end
    end

standard_deviation(array, &block)
Alias for: std
std(array, &block) click to toggle source

Standard deviation of a sample.

    # File lib/facets/math/std.rb, line 7
    def self.std(array, &block)
      sqrt(variance(array, &block))
    end

Also aliased as: standard_deviation
stderr(array) click to toggle source

Calculates the standard error of a sample.

    # File lib/facets/math/std.rb, line 22
    def self.stderr(array)
      return 0.0 if array.size < 2
      std(array) / sqrt(array.size)
    end

sum(array) { || ... } click to toggle source

Returns sum. When a block is given, summation is taken over the each result of block evaluation.

    # File lib/facets/math/sum.rb, line 6
    def self.sum(array) #:yield:
      sum = 0.0
      if block_given?
        array.each{|i| sum += yield(i)}
      else
        array.each{|i| sum += i}
      end
      sum
    end

summed_sqdevs(array) click to toggle source

The sum of the squared deviations from the mean.

    # File lib/facets/math/summed_sqdevs.rb, line 8
    def self.summed_sqdevs(array)
      return 0 if array.size < 2
      m = mean(array)
      sum(array.map{ |x| (x - m) ** 2 })
    end

tgamma(x) click to toggle source

Exp of LGamma.

    # File lib/facets/math/tgamma.rb, line 6
    def self.tgamma(x)
      exp(lgamma(x))
    end

theil_index(array) click to toggle source

Calculates the Theil index (a statistic used to measure economic inequality).

TI = sum_{i=1}^N frac{x_i}{sum_{j=1}^N x_j} ln frac{x_i}{bar{x}}

    http://en.wikipedia.org/wiki/Theil_index

    # File lib/facets/math/theil_index.rb, line 14
    def self.theil_index(array)
      return -1 if array.size <= 0 or any? { |x| x < 0 }
      return  0 if array.size <  2 or all? { |x| approx_equal(x, 0) }
      m = mean(array)
      s = sum(array).to_f
      inject(0) do |theil, xi|
              theil + ((xi > 0) ? (log(xi.to_f/m) * xi.to_f/s) : 0.0)
      end
    end

variance(array, &block) click to toggle source

    # File lib/facets/math/variance.rb, line 6
    def self.variance(array, &block)
      sum2 = if block_given?
        sum(array){ |i| j = block[i]; j*j }
      else
        sum(array){ |i| i**2 }
      end
      sum2/array.size - mean(array, &block)**2
    end

variance2(array) click to toggle source

Variance of the sample. Variance of 0 or 1 elements is 0.0.

TODO: Same as variance? Then choose one.

    # File lib/facets/math/variance.rb, line 19
    def self.variance2(array)
      return 0.0 if array.size < 2
      summed_sqdevs(array) / (array.size - 1)
    end

### Public Instance Methods

acsch(x) click to toggle source

Area cosecans hyperbolicus of `x`

    # File lib/facets/math/acsch.rb, line 4
    def acsch(x)
      ::Math.log(1.0 / x + Math.sqrt(1.0 + 1.0 / (x * x)))
    end

asech(x) click to toggle source

Area secans hyperbolicus of `x`

    # File lib/facets/math/asech.rb, line 4
    def asech(x)
      log((1.0 + sqrt(1.0 - x * x)) / x)
    end