A little diversion (probability fun in Ruby)

[Daniel Martin]

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.chooseboy, :girl)
child2 = u.chooseboy, :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.chooseboy, :girl), u.chooseboy, :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.sendsucceed)
else
return ptree.sendnormalize)
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.chooseboy, :girl)
child2 = u.chooseboy, :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

 

[Matt Todd]

This sounds really cool. I regret not paying more attention in my
statistics class. I've been rereading my textooks on math, and soon
enough I'll actually understand how to produce probability! Hah.

Thanks for this cool code. I'll probably play with it while I learn
about probability.

Cheers, and good luck with that interview!

M.T.