Class: AliasTable

Inherits:
Object
  • Object
show all
Defined in:
lib/aliastable.rb

Overview

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.

Instance Method Summary collapse

Constructor Details

#initialize(x_set, p_value) ⇒ AliasTable

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_values are different lengths.

  • RuntimeError if any p_value is negative.

  • RuntimeError if p_values don’t sum to one. Rationals will avoid this.



28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
 
# File 'lib/aliastable.rb', line 28

def initialize(x_set, p_value)
  if x_set.length != p_value.length
    fail 'Args to AliasTable must be vectors of the same length.'
  end
  fail 'p_values must be positive' unless p_value.all? { |value| value > 0 }
  @p_primary = p_value.map(&:to_r)
  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

Instance Method Details

#generate ⇒ Object

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.



59
60
61
62
 
# File 'lib/aliastable.rb', line 59

def generate
  column = rand(@x.length)
  rand <= @p_primary[column] ? @x[column] : @alias[column]
end