class AliasTable

Generate values from a categorical distribution in constant time, regardless of the number of categories. This clever algorithm uses conditional probability to construct a table comprised of columns which have a primary value and an alias. Generating a value consists of picking any column (with equal probabilities), and then picking between the primary and the alias based on appropriate conditional probabilities.

Public Class Methods

new(x_set, p_value) click to toggle source

Construct an alias table from a set of values and their associated probabilities. Values and their probabilities must be synchronized, i.e., they must be arrays of the same length. Values can be anything, but the probabilities must be positive Rational numbers that sum to one.

Arguments

  • x_set -> the set of values from which to generate.

  • p_value -> the synchronized set of probabilities associated with the value set. These values should be Rationals to avoid rounding errors.

Raises

  • RuntimeError if x_set and +p_value+s are different lengths.

  • RuntimeError if any p_value is negative.

  • RuntimeError if +p_value+s don't sum to one. Rationals will avoid this.

# File lib/aliastable.rb, line 28
def initialize(x_set, p_value)
  fail 'x_set & p_value must have same length.' if x_set.size != p_value.size
  fail 'p_values must be positive' unless p_value.all? { |value| value > 0 }
  @p_primary = p_value.map(&:rationalize)
  fail 'p_values must sum to 1' unless @p_primary.reduce(:+) == Rational(1)
  @x = x_set.clone.freeze
  @alias = Array.new(@x.length)
  parity = Rational(1, @x.length)
  group = @p_primary.each_index.group_by { |i| @p_primary[i] <=> parity }
  deficit_set = group[-1]
  surplus_set = group[1]
  until deficit_set.empty?
    deficit = deficit_set.pop
    surplus = surplus_set.pop
    @p_primary[surplus] -= parity - @p_primary[deficit]
    @p_primary[deficit] /= parity
    @alias[deficit] = @x[surplus]
    if @p_primary[surplus] == parity
      @p_primary[surplus] = Rational(1)
    else
      (@p_primary[surplus] < parity ? deficit_set : surplus_set) << surplus
    end
  end
end

Public Instance Methods

generate() click to toggle source

Return a random outcome from this object's distribution. The generate method is O(1) time, but is not an inversion since two uniforms are used for each value that gets generated.

# File lib/aliastable.rb, line 57
def generate
  column = rand(@x.length)
  rand <= @p_primary[column] ? @x[column] : @alias[column]
end