From: Daniel Martin <martin@...>
Date: 2006-08-05T00:27:04+09:00
Subject: A little diversion (probability fun in Ruby)

So I can't sleep this morning, and I've got a big interview with
$PRESTIGOUS_COMPANY that I need to get my mind off of, so I thought
I'd share a little thing I did in Ruby a few months ago with the list.

If you've never heard of the "Monty Hall" problem, suffice it to say
that it's a big famous problem in probability
(http://en.wikipedia.org/wiki/Monty_Hall_problem) that is sometimes
used to see if people really, really know their probability.  Anyway,
I've written a ruby class that allows one to investigate probabilities
such as the monty hall problem.  For example:

puts('===Monty Hall, classic version===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  guessdoor = u.choose(1,2,3)
  remaining_doors = [1,2,3].select{ |x|
     x != treasuredoor and x != guessdoor }
  showdoor = u.choose(*remaining_doors)
  if (treasuredoor == guessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

Produces:

===Monty Hall, classic version===
You should switch
==>     2/3
You should stay
==>     1/3

Which is in fact the known probabilities.  While playing with this, I
also have discovered subtleties to the problem, such as this
variation:

# variation where the host doesn't know where the treasure is either,
# and sometimes accidentally reveals the treasure.  In that case,
# the producers declare a do-over, and the segment is retaped.
puts('===Monty Hall with ignorant host and do-overs===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  guessdoor = u.choose(1,2,3)
  remaining_doors = [1,2,3].select{|x| x != guessdoor}
  showdoor = u.choose(*remaining_doors)
  u.assert(treasuredoor != showdoor)  # i.e. we're not a do-over
  if (treasuredoor == guessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

Produces:

===Monty Hall with ignorant host and do-overs===
You should switch
==>     1/2
You should stay
==>     1/2

This is in fact also known to be the correct result.

I can tackle other probability puzzles/games, such as:

# A family has two children.  One of them is a boy.
# What is a chance that the other is a boy?
puts "===Boy-Girl problem==="
ProbabilityTree.runreport(1.to_r) { |u|
  child1 = u.choose(:boy, :girl)
  child2 = u.choose(:boy, :girl)
  u.assert(child1 == :boy || child2 == :boy)
  if child1 == :boy and child2 == :boy then
    u.report "Both children are boys"
  else
    u.report "Other child is a girl"
  end
}.display

===Boy-Girl problem===
Other child is a girl
==>     2/3
Both children are boys
==>     1/3

Anyway, I've found this fun to play around with, and useful in
clarifying many subtleties in probability problems and how to work
them.  It takes a bit of practice to translate properly from problem
statement to code, and along the way gets you to think about exactly
what's going on.

For example, the boy-girl problem could be incorrectly coded as:

ProbabilityTree.runreport(1.to_r) { |u|
  children = [u.choose(:boy, :girl), u.choose(:boy, :girl)]
  choose_index = u.choose(0,1)
  u.assert(children[choose_index] == :boy)
  u.report "Other child is a #{children[1-choose_index]}"
}.display

Which is a different problem.  (This problem is: "A couple has two
children, and you pick one at random.  That child is a boy.  What is
the chance that the other child is also a boy?")

You can also choose things with unequal probability, so you can
simulate loaded dice if you want to.  I tried to simulate the
probabilities of blood types for my offspring given everything I know
about blood type distributions in the US and among my relatives, but
that never got anywhere in ruby; it's just a bit too slow.  I do have
this same idea implemented in Haskell, however, and that was able to
produce an answer to the blood-type question in a reasonable amount of
time.

Anyway, here's the code:

____________probworlds.rb_________________
# whatever
class Amb
  class ExhaustedError < RuntimeError; end

  def initialize
    @fail = [proc { fail ExhaustedError, "amb tree exhausted" }]
  end

  def choose(*choices)
    enum_choose(choices)
  end
  
  def enum_choose(choices)
    choices.each { |choice|
      callcc { |fk|
        @fail << fk
        if choice.respond_to? :call
          return choice.call
        else
          return choice
        end
      }
    }
    @fail.pop.call
  end

  def failure
    choose
  end

  def assert(cond)
    failure unless cond
  end
end

class ProbabilityTree
  private_class_method :new

  private

  def initialize(amb,one=nil)
    @amb = amb
    @reporthash = Hash.new(0)
    @currentreport = nil
    @currentprob = (one || 1)
    @one = one || 1.0
    @probtotal = 0
  end
    
  def save
    [ @currentreport, @currentprob ]
  end
    
  def restore(savedstate)
    @currentreport, @currentprob = savedstate
  end
    
  def normalize
    @reporthash.each_key { |k| @reporthash[k] /= @probtotal }
    @probtotal = 1 unless @reporthash.empty?
    @reporthash
  end
    
  def succeed
    @reporthash[@currentreport] += @currentprob
    @probtotal += @currentprob
    @amb.failure
  end
    
  public
    
  def failure
    @amb.failure
  end
    
  def choose_h(h)
    choose_p(*(h.to_a))
  end
    
  def choose(*rest)
    choose_p(*(rest.map {|x| [x, @one/rest.size]}))
  end
    
  def choose_p(*args)
    savedstate = save()
    x, prob = @amb.enum_choose(args)
    restore(savedstate)
    @currentprob *= prob
    return x
  end
    
  def report(string)
    if (@currentreport and @currentreport != '')
      @currentreport = @currentreport.to_s
      @currentreport += "\n#{string.to_s}"
    else
      @currentreport = string
    end
  end
    
  def assert(cond)
    failure unless cond
  end
    
  def ProbabilityTree.universe(one=nil)
    ae = Amb.new
    ptree = new(ae,one)
    if (ae.choose(true,false))
      yield ptree
      ptree.send(:succeed)
    else
      return ptree.send(:normalize)
    end
  end

  def ProbabilityTree.report(h)
    rval = ""
    h.keys.sort {|x,y| h[y] <=> h[x]}.each { |k|
      rval += "#{k}\n==>\t#{h[k]}\n"
    }
    rval
  end

  def ProbabilityTree.runreport(one=nil)
    report(universe(one) {|u| yield u})
  end
end

# reference http://en.wikipedia.org/wiki/Three_cards_problem
puts 'Three-card problem:'
ProbabilityTree.runreport { |u|
  card = u.choose('bw', 'ww', 'bb')
  side = u.choose(0,1)
  u.assert(card[side] == ?b);
  u.report "Other side is #{card[1-side,1]}"
}.display

puts 'Three-card problem with someone who always shows you the black
side if possible:'
ProbabilityTree.runreport { |u|
  card = u.choose('bw', 'ww', 'bb')
  side = u.choose(0,1)
  if card[side] != ?b and card[1-side] == ?b then
    side = 1-side
  end
  u.assert(card[side] == ?b);
  u.report "Other side is #{card[1-side,1]}"
}.display

# Monty Hall
require "rational"

# standard version
puts('===Monty Hall, classic version===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  guessdoor = u.choose(1,2,3)
  remaining_doors = [1,2,3].select{|x|
    x != treasuredoor and x != guessdoor }
  showdoor = u.choose(*remaining_doors)
  if (treasuredoor == guessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

# variation where the host doesn't know where the treasure is either,
# and sometimes accidentally reveals the treasure.  In that case,
# the producers declare a do-over, and the segment is retaped.
puts('===Monty Hall with ignorant host and do-overs===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  guessdoor = u.choose(1,2,3)
  remaining_doors = [1,2,3].select{|x| x != guessdoor}
  showdoor = u.choose(*remaining_doors)
  u.assert(treasuredoor != showdoor)  # i.e. we're not a do-over
  if (treasuredoor == guessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

# The lesson?
# use assert for stuff that "just happened".  Pass a
# different list to choose for stuff that was deliberate.

# You and Bob both choose doors.  Monty Hall tells one of you
# to go home, and opens the losers door, because they chose poorly.
# Should the remaining person switch doors or stick with their
# original guess?
puts('===Monty Hall with Bob===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  myguessdoor = u.choose(1,2,3)
  bobguessdoor = u.choose(*([1,2,3].select{|i| i != myguessdoor}))
  elimchoices = [[:me,myguessdoor],[:bob,bobguessdoor]]
  elimchoices = elimchoices.select{|i| i[1] != treasuredoor}
  elimchoice = u.choose(*elimchoices)
  u.assert(elimchoice[0] == :bob)
  if (treasuredoor == myguessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

# Same as above, but the host hates Bob and will eliminate
# him early whenever possible.  Now, if Bob is eliminated
# what should you do?
puts('===Monty Hall hates Bob===')
ProbabilityTree.runreport(1.to_r) { |u|
  treasuredoor = u.choose(1,2,3)
  myguessdoor = u.choose(1,2,3)
  bobguessdoor = u.choose(*([1,2,3].select{|i| i != myguessdoor}))
  elimchoices = [[:me,myguessdoor],[:bob,bobguessdoor]]
  elimchoices = elimchoices.select{|i| i[1] != treasuredoor}
  elimchoice = nil
  if elimchoices.size == 2 then
    # If we're both still in the running to be eliminated, choose Bob
    elimchoice = elimchoices[1]
  else
    # Choose whoever is left
    elimchoice = elimchoices[0]
  end
  u.assert(elimchoice[0] == :bob)
  if (treasuredoor == myguessdoor)
    u.report "You should stay"
  else
    u.report "You should switch"
  end
}.display

# This doesn't work
puts "Sum - naive"
ProbabilityTree.runreport(1.to_r) { |u|
  sum = 0
  3.times {
    i = u.choose(0,1)
    sum += i
  }
  u.report sum
}.display

# This does work
puts "Sum - store state"
ProbabilityTree.runreport(1.to_r) { |u|
  sum = 0
  3.times {
    psum = sum
    i = u.choose(0,1)
    sum = psum + i
  }
  u.report sum
}.display

# This works too
puts "Sum - hidden store state"
ProbabilityTree.runreport(1.to_r) { |u|
  sum = 0
  3.times {
    sum += u.choose(1,0)
  }
  u.report sum
}.display

# and this
puts "Sum - choose state"
ProbabilityTree.runreport(1.to_r) { |u|
  sum = 0
  3.times {
    i,sum = u.choose([0,sum],[1,sum])
    sum += i
  }
  u.report sum
}.display

# and this
puts "Sum - inject"
ProbabilityTree.runreport(1.to_r) { |u|
  sum = (0 ... 3).inject(0) { |s,|
    s + u.choose(0,1)
  }
  u.report sum
}.display

# A family has two children.  One of them is a boy.
# What is a chance that the other is a boy?
puts "===Boy-Girl problem==="
ProbabilityTree.runreport(1.to_r) { |u|
  child1 = u.choose(:boy, :girl)
  child2 = u.choose(:boy, :girl)
  u.assert(child1 == :boy || child2 == :boy)
  if child1 == :boy and child2 == :boy then
    u.report "Both children are boys"
  else
    u.report "Other child is a girl"
  end
}.display