Module: ChrisMath

Includes:
Math
Defined in:
lib/chris_lib/chris_math.rb

Overview

Numerically focused helpers and distribution utilities.

Instance Method Summary collapse

Instance Method Details

#bi_gaussian_rand ⇒ Array<Float>

Generates a pair of independent standard normal deviates

Returns:

  • (Array<Float>) —

    two-element array [z0, z1]



234
235
236
237
238
239
240
 
# File 'lib/chris_lib/chris_math.rb', line 234

def bi_gaussian_rand
  u1 = rand()
  u2 = rand()
  z0 = sqrt(-2*log(u1))*cos(2*PI*u2)
  z1 = sqrt(-2*log(u1))*sin(2*PI*u2)
  [z0,z1]
end

#combinatorial(n, k) ⇒ Integer

Binomial coefficient “n choose k”

Parameters:

  • n (Integer)
  • k (Integer)

Returns:

  • (Integer)

Raises:

  • (ArgumentError) 


319
320
321
322
323
324
325
 
# File 'lib/chris_lib/chris_math.rb', line 319

def combinatorial(n, k)
  #from rosetta code
  raise ArgumentError, 'n must be a non-negative Integer' unless n.is_a?(Integer) && n >= 0
  raise ArgumentError, 'k must be a non-negative Integer' unless k.is_a?(Integer) && k >= 0
  raise ArgumentError, 'k must be <= n' if k > n
  (0...k).inject(1) do |m,i| (m * (n - i)) / (i + 1) end
end

#combinatorial_distribution(n, r, p) ⇒ Float

Cumulative probability of at least r successes in n Bernoulli trials

Parameters:

  • n (Integer)
  • r (Integer)
  • p (Float) —

    probability of success per trial

Returns:

  • (Float)

Raises:

  • (RuntimeError) —

    when r is outside 0..n



284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
 
# File 'lib/chris_lib/chris_math.rb', line 284

def combinatorial_distribution(n,r,p)
  raise ArgumentError, 'n must be a non-negative Integer' unless n.is_a?(Integer) && n >= 0
  raise ArgumentError, 'r must be a non-negative Integer' unless r.is_a?(Integer) && r >= 0
  raise ArgumentError, 'p must be between 0 and 1' unless p.is_a?(Numeric) && p.between?(0.0,1.0)
  # probability that r out n or greater hits with the 
  # probability of one hit is p.
  if r <= n && r >= 0
  sum = 0
  (r..n).each do |k|
    sum += combinatorial(n,k)*(p)**(k)*(1-p)**(n - k)
  end
  sum
  else
    raise 'Error, #{r} must be >= 0  and <= #{n}'
  end
end

#factorial(n) ⇒ Integer

Note:

Warns and returns ‘nil` when `n` is negative.

Factorial without overflow protection

Parameters:

  • n (Integer)

Returns:

  • (Integer)

Raises:

  • (ArgumentError) 


305
306
307
308
309
310
311
312
313
 
# File 'lib/chris_lib/chris_math.rb', line 305

def factorial(n)
  #from rosetta code
  raise ArgumentError, 'n must be an Integer' unless n.is_a?(Integer)
  if n.negative?
    warn 'factorial called with negative n; returning nil'
    return nil
  end
  (1..n).inject {|prod, i| prod * i}
end

#gaussian_array(n = 1) ⇒ Array<Float>

Note:

Returns an empty array and warns when ‘n <= 0`.

Generates normally distributed random numbers using Box-Muller

Parameters:

  • n (Integer) (defaults to: 1) —

    number of samples to generate

Returns:

  • (Array<Float>) —

    array of N(0,1) samples

Raises:

  • (ArgumentError) 


219
220
221
222
223
224
225
226
227
228
229
230
 
# File 'lib/chris_lib/chris_math.rb', line 219

def gaussian_array(n = 1)
  raise ArgumentError, 'n must be an Integer' unless n.is_a?(Integer)
  if n <= 0
    warn 'gaussian_array called with non-positive sample size; returning empty array'
    return []
  end
  (1..n).map do
    u1 = rand()
    u2 = rand()
    sqrt(-2*log(u1))*cos(2*PI*u2)
  end
end

#gaussian_rand(mean, std) ⇒ Float

Note:

Warns when ‘std <= 0`, returning the mean or using the absolute value.

Generates a normally distributed random number with the provided mean and std

Parameters:

  • mean (Numeric)
  • std (Numeric)

Returns:

  • (Float)

Raises:

  • (ArgumentError) 


247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
 
# File 'lib/chris_lib/chris_math.rb', line 247

def gaussian_rand(mean,std)
  raise ArgumentError, 'mean must be Numeric' unless mean.is_a?(Numeric)
  raise ArgumentError, 'std must be Numeric' unless std.is_a?(Numeric)
  if std.negative?
    warn 'std supplied to gaussian_rand was negative; using absolute value'
    std = std.abs
  elsif std.zero?
    warn 'std supplied to gaussian_rand was zero; returning mean'
    return mean
  end
  u1 = rand()
  u2 = rand()
  z0 = sqrt(-2*log(u1))*cos(2*PI*u2)
  z0*std + mean
end

#std(values) ⇒ Float

Sample standard deviation for raw values

Parameters:

  • values (Array<Numeric>)

Returns:

  • (Float)

Raises:

  • (RuntimeError) —

    when fewer than two observations are provided



268
269
270
271
272
273
274
275
 
# File 'lib/chris_lib/chris_math.rb', line 268

def std(values)
  raise ArgumentError, 'values must respond to #length and #inject' unless values.respond_to?(:length) && values.respond_to?(:inject)
  n = values.length
  raise 'n = #{n} but must be greater than 1' if n < 2
  m = mean(values)
  sum = values.inject { |s,v| s + (v**2 - m**2)}
  sqrt(sum.to_f/(n - 1))
end