class RV::Gamma

Gamma generator based on Marsaglia and Tsang method Algorithm 4.33

Produces gamma RVs with expected value alpha * beta.

Arguments
  • alpha -> the shape parameter (alpha > 0; default: 1).

  • beta -> the rate parameter (beta > 0; default: 1).

  • rng -> the (Enumerable) source of U(0, 1)'s (default: U_GENERATOR)

Attributes

alpha[R]
beta[R]

Public Class Methods

new(alpha: 1.0, beta: 1.0, rng: U_GENERATOR) click to toggle source
# File lib/random_variates.rb, line 280
def initialize(alpha: 1.0, beta: 1.0, rng: U_GENERATOR)
  raise 'Alpha and beta must be positive.' if alpha <= 0 || beta <= 0

  @alpha = alpha
  @beta = beta
  @rng = rng
  @std_normal = Normal.new(rng: rng)
end

Public Instance Methods

next() click to toggle source
# File lib/random_variates.rb, line 289
def next
  __gen__(@alpha, @beta)
end

Private Instance Methods

__gen__(alpha, beta) click to toggle source
# File lib/random_variates.rb, line 295
def __gen__(alpha, beta)
  if alpha > 1
    z = v = 0.0
    d = alpha - 1.0 / 3.0
    c = (1.0 / 3.0) / Math.sqrt(d)
    loop do
      loop do
        z = @std_normal.next
        v = 1.0 + c * z
        break if v > 0
      end
      z2 = z * z
      v = v * v * v
      u = @rng.next
      break if u < 1.0 - 0.0331 * z2 * z2
      break if Math.log(u) < (0.5 * z2 + d * (1.0 - v + Math.log(v)))
    end
    d * v * beta
  else
    __gen__(alpha + 1.0, beta) * (@rng.next**(1.0 / alpha))
  end
end