Likelihood weighting and sampling

sampling.rb

# sampling.rb -- approximate the probability distribution of Pr[A,B|-C,+F]
#
# Author: Jon-Michael Deldin

def label_table(label, n, fields)
  puts "#{label} (n = #{n}):"
  hdr = "%-15s %-10s" % fields
  puts hdr
  puts "-" * hdr.length
end

def print_row(sample, prob)
  puts "%-15s %0.3f" % [sample, prob]
end

# Switch "+A" to "-A"
def swap_sign(node)
  sign = node.chars.to_a[0] == '+' ? '-' : '+'
  sign + node[1..-1]
end

# Given probabilities
Pr = {
  '+A'      => 0.1,
  '+B'      => 0.7,
  '+C|+A'   => 0.3,
  '+C|-A'   => 0.8,
  '+D|+A+B' => 0.3,
  '+D|-A+B' => 0.8,
  '+D|+A-B' => 0.1,
  '+D|-A-B' => 0.2,
  '+F|+D'   => 0.9,
  '+F|-D'   => 0.8,
}

# Generate complementary probabilities
Pr.keys.each do |k|
  Pr[swap_sign(k)] = 1 - Pr[k]
end

# Select =node= or -node depending on the Pr[node]
def Pr.sample_node(node)
  rand <= self[node] ? node : swap_sign(node)
end

# Same as sample_node, but handles given clauses
def Pr.sample_node_with_given(node)
  chr = node.scan(/^[+-][A-Z]/).first
  rand <= Pr[node] ? chr : swap_sign(chr)
end

# Returns a hash of various AB combinations we're interested in
def generate_sum_hash
  {'+A+B' => 0, '+A-B' => 0, '-A+B' => 0, '-A-B' => 0}
end

def likelihood_weighting(n)
  sums  = generate_sum_hash
  total = 0.0

  n.times do
    a = Pr.sample_node '+A'
    b = Pr.sample_node '+B'
    d = Pr.sample_node_with_given "+D|#{a}#{b}"

    total_weight = Pr["-C|#{a}"] * Pr["+F|#{d}"]
    sums["#{a}#{b}"] += total_weight
    total += total_weight
  end

  label_table('Likelihood Weighting', n, %w(Sample Weight))
  sums.each do |pair, weight|
    print_row("Pr[#{pair}|-C+F]", weight/total)
  end
end

def rejection_sampling(n)
  kept  = generate_sum_hash
  total = 0

  n.times do
    a = Pr.sample_node '+A'
    b = Pr.sample_node '+B'
    c = Pr.sample_node_with_given "+C|#{a}"
    d = Pr.sample_node_with_given "+D|#{a}#{b}"
    f = Pr.sample_node_with_given "+F|#{d}"

    next if c[0] == '-' || f[0] == '+'

    key = "#{a}#{b}"
    kept[key] += 1
    total     += 1
  end

  label_table('Rejection Sampling', "#{n}, actual = #{total}", %w(Sample Prob))
  kept.each do |key, kept_ct|
    print_row("Pr[#{key}|-C,+F]", kept_ct.to_f/total)
  end
end

if $0 == __FILE__
  n = 10_000
  likelihood_weighting(n)
  puts
  rejection_sampling(n)
end

% ruby sampling.rb
Likelihood Weighting (n = 10000):
Sample          Weight    
--------------------------
Pr[+A+B|-C+F]   0.190
Pr[+A-B|-C+F]   0.074
Pr[-A+B|-C+F]   0.535
Pr[-A-B|-C+F]   0.201

Rejection Sampling (n = 10000, actual = 1024):
Sample          Prob      
--------------------------
Pr[+A+B|-C,+F]  0.035
Pr[+A-B|-C,+F]  0.019
Pr[-A+B|-C,+F]  0.585
Pr[-A-B|-C,+F]  0.361