Monte Carlo Solutions to 50 Challenging Problems in Probabilities Blog
When I learned about Mosteller's Fifty Challenging Problems in Probabilities, I had to have a copy. I'm a big fan of so-called "recreational mathematics," and was working on honing my statistical intuitions at the time, so it was right up my alley. When I dug into it, though, I found that the statistics were… involved. Looking at the problems, it was clear that most of them were trivially solvable in simulations. This page is the outcome of that.
The biggest takeaway from my experience here is this: when you face a probabilistic problem, you should start by finding a way to simulate the behavior. Only reach for statistics when you're checking your work on that simulation. There are some gotchas out there: for instance, you need to know what kind of probability distribution your model should use (ie is it normal, uniform, or something esoteric?). However, a bit of experimentation will usually suss these things out.
[As I prepare this for (re-)publication, I've "grown up" a bit in my statistical sophistication. I still believe in the above, but would likely reach for different tools today. In particular, numpy with ipython and GNU R are both powerful environments for the kinds of exploratory work done below. The irony is that, had I used them, I probably wouldn't have these kind of concrete examples to share with you all. –jbm 20130127 ]
Problem 2 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 2: To encourage Elmer's promising tennis carrer, his father offers him a prize if he wins two or more sets in a row in a three-set series to be played with this father and the club champion alternately: father-champion-father or champ-dad-champ. The champion is a better player than Elmer's father. Which order should Elmer choose?
#!/usr/bin/env ruby
# This is another search problem. We're not entirely
# sure if the relative gap between the pro and the
# father matters, and we're not even sure that the son
# can beat the father, much less the pro. Alternately,
# he might beat them both handily.
#
# So, we'll take sample points: we'll give the son 11
# different chances against the pro: 0% chance of
# winning, 10% chance of winning, ... , 90% chance,
# 100%. Then, for each of these, we'll give him all
# the same chances up to this chance against his father.
# At each of these, we run a number of trials with both
# orders.
#
# Yes, this gets expensive: it's the upper triangle of a
# 11x11 matrix, so we're going to do
# 11x11/2*2*TRIALS = 121*TRIALS trials.
# But, hey, cycles are cheap.
TRIALS=10000
RANGE=[0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
n_dpd = 0
n_pdp = 0
def win(series)
did_win = 0 # 1 is "won last", 2 is "won bet"
series.each do |p_now|
if p_now > rand()
did_win += 1
return true if did_win == 2
else
did_win = 0
end
end
return false
end
RANGE.each do |p_pro|
RANGE.each do |p_dad|
break if p_dad > p_pro # we'll let dad == pro through for now
TRIALS.times do
n_dpd += 1 if (win([p_dad, p_pro, p_dad]))
n_pdp += 1 if (win([p_pro, p_dad, p_pro]))
end
end
end
puts "Total Wins: dad/pro/dad=#{n_dpd}, pro/dad/pro=#{n_pdp}"
Problem 3 of Monte Carlo solutions to Fifty Challenging Problems...
In a jury, each individual member has probability p of making the right decision. Which is more likely to find the correct decision: a 3-person jury in which the third member decides by flipping a coin (50/50), or a 1-person jury? (NB: these juries are majority rule, not unanimous)
#!/usr/bin/env ruby
# First off, we need to model the juries:
# 1. One-man jury: just p every time
# 2. Three-man jury: p,p,0.5 (majority rules)
#
# The goal is: which jury is more likely to be correct?
# Apparent sexism note: if "woman" was shorter than "man",
# these would all be n_woman_judgement. I'm lazy.
def one_man_judgement(p)
sample = rand()
if (sample < p)
return true
else
return false
end
end
def three_man_judgement(p)
j1 = one_man_judgement(p)
j2 = one_man_judgement(p)
j3 = one_man_judgement(0.5)
if (j1 && j2 || j1 && j3 || j2 && j3)
return true
else
return false
end
end
# So, now we can get the judgement of a jury by calling one of the above.
# if it's true, they got the correct judgement; if false, they got it wrong.
# Since we're leaving p unspecified, we should sample it at multiple
# values. I usually go for 20 right up front, since it's often over
# the precision needs of a problem. If you see anything fishy (esp
# non-monotonic behavior), you probably want to run again with tighter
# buckets, to get a better idea of the shape of the curve.
SAMPLE_WIDTH=0.05
TRIALS=10000
one_man_wins = 0
three_man_wins = 0
p = 0.0
while (p <= 1.0)
# now, we create the confusion matrix:
one_man_outcomes = { :correct => 0, :incorrect => 0 }
three_man_outcomes = { :correct => 0, :incorrect => 0 }
TRIALS.times do
o = one_man_judgement(p)
t = three_man_judgement(p)
one_man_outcomes[ o ? :correct : :incorrect ] += 1
three_man_outcomes[ t ? :correct : :incorrect ] += 1
end
p_one = one_man_outcomes[:correct] / TRIALS.to_f
p_three = three_man_outcomes[:correct] / TRIALS.to_f
puts "==== #{p} // #{TRIALS} ===="
puts "sample\ttrue\tfalse\tP"
puts "one\t#{one_man_outcomes[:correct]}\t#{one_man_outcomes[:incorrect]}\t#{p_one}"
puts "three\t#{three_man_outcomes[:correct]}\t#{three_man_outcomes[:incorrect]}\t#{p_three}"
one_man_wins += one_man_outcomes[:correct]
three_man_wins += three_man_outcomes[:correct]
p += SAMPLE_WIDTH
end
puts
puts "FINAL: At #{one_man_wins} samples one-man is better; at #{three_man_wins} three-man is better"
Problem 5 of Monte Carlo solutions to Fifty Challenging Problems...
Draw a grid of 1" squares on a table, then toss a penny (diameter 3/4") onto it. What's the probability the penny doesn't cross any line on the table? (Now, turn this into a carnival game and make some money.)
#!/usr/bin/env ruby
TRIALS=10000
# The goal in this one is to determine the probability
# that the coin is inside the square. We can model
# this one of two ways: the probability of crossing
# a line, or the probability of not crossing any lines.
#
# P(crossing a line) is pretty hard, with four different
# cases to consider (unless you model the grid in a
# slightly inside-out way. After going through the below,
# think of how you could rework the setup to model crossing
# either line in the same manner.)
#
# Instead, we model the problem as
# P(not crossing a line)
#
# We model the coin as a point (x,y), with a circle
# drawn 3/8" around it, yielding a diameter of 3/4".
#
# Thus, to win, you must have x and y in the
# range [3/8,5/8].
#
RANGE_LOW = 3.0/8.0
RANGE_HIGH = 5.0/8.0
def win()
x = rand()
y = rand()
return ( (x > RANGE_LOW) && (x < RANGE_HIGH) &&
(y > RANGE_LOW) && (y < RANGE_HIGH) )
end
wins = 0
TRIALS.times {
wins +=1 if win()
}
puts "Out of #{TRIALS} coin tosses, we won #{wins}."
puts "The probability of winning is approximately #{wins.to_f/TRIALS.to_f}"
Problem 6 of Monte Carlo solutions to Fifty Challenging Problems...
What is the expected loss per unit stake in Chuck-a-Luck? Only single-die bets are allowed.
#!/usr/bin/env ruby
TRIALS=50000
payout = 0
# First off, the game rules. If we get one die right,
# it pays out 1 times our stake. If we get two dice,
# twice the stake, and three dice yields three times
# the stake. This function returns the multiplicand
# of our stake; if we lose, it returns 0.
#
def win(bet)
ret = 0
d1 = 1+rand(6)
d2 = 1+rand(6)
d3 = 1+rand(6)
ret += 1 if (d1 == bet)
ret += 1 if (d2 == bet)
ret += 1 if (d3 == bet)
ret +=1 if ret > 0 # we get our stake back, too
return ret
end
# So, now, we play the game. Since the dice are fair,
# all numbers are equally likely, and we can bet on any
# number we choose each game. To simplify things,
# we'll always bet on 3.
TRIALS.times do
payout += win(3) # we bet 1 unit on rolling a three
end
t = TRIALS # make the prints fit.
puts "After #{t} rounds of Chuck-a-Luck, we spent $#{t}"
puts "We got back $#{payout}. Thus, our expected loss"
puts " per unit stake is $#{(TRIALS-payout)/TRIALS.to_f}."
Problem 7 of Monte Carlo solutions to Fifty Challenging Problems...
Mr Brown has a problem: he compulsively puts a dollar on number 13 in American roulette. To help him see the error of his ways, Mr Pink gives him a bet: Pink bets Brown $20 if, after 36 rounds of roulette, Brown is ahead. (That is: if Brown is up after 36 compulsive rounds, he gets an extra $20; if he's behind, he loses an extra $20.)
#!/usr/bin/env ruby
# We'll play from Mr Brown's perspective.
#
# Modelling this game is tricky. We need to play 36
# rounds of roulette, then see if we're up. If so,
# we get an extra $20; if not, we lose $20. At the
# end, we want the expected payout per unit stake.
TRIALS=10000
def payout()
p = 0
36.times do
p += 36 if ( rand() < 1/38.to_f)
end
if p > 26
p += 20
else
p -= 20
end
# puts "This round: $36 in, $#{p} out\n"
p
end
total_payout = 0
total_bet = 0
TRIALS.times {
total_payout += payout()
total_bet += 36
}
puts "After #{TRIALS} rounds of betting, we bet #{total_bet}"
puts "We won #{total_payout}, so the expected loss is "
puts " $#{(total_bet-total_payout)/total_bet.to_f}"
Problem 8 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 8: What is the chance of being dealt a perfect hand in bridge? (being dealt 13 cards of the same suit?)
This is one place where Monte Carlo falls flat on its face. Modeling this in any sort of reasonably simple (read: naive) fashion results in a program that takes quite some time to run. It's far better to find the probability through direct combinatorial calculations. This still poses some interesting programming problems, but they're of the numerical sort, so they won't appear in this series. (It would be really great if someone would post a comment linking to a solution of this problem using hadoop on EC2 or EMR!)
#!/usr/bin/env ruby
TRIALS=100
# We model this by laying out the 52 cards in an
# array, calling the spades 1..13.
#
# Then, we "deal" by breaking it up into contiguous
# 13-card chunks. If any the 4 chunks of 13 are
# 1..13, we win the "deal a perfect hand" game.
# Otherwise, we lose.
perfect_hands = 0
# I had some concerns about the performance here.
def perfect_hand_by_sort?(a)
ap = a.sort
return 12 == ap[12]-ap[0]
end
def perfect_hand_by_scan?(a, handno)
x = (handno-1)*13
min = 0
max = 53
while (x < (handno)*13)
ai = a[x]
min = ai if ai < min
max = ai if ai > max
x += 1
end
return 12 == max-min
end
def perfect_hand_by_minmax?(a, handno)
x = (handno-1)*13
min = a[x..x+12].min
max = a[x..x+12].max
return 12 == max-min
end
t_by_sort = 0.0
t_by_scan = 0.0
t_by_minmax = 0.0
TRIALS.times do
cards = (1..52).to_a
cards.shuffle!
t0 = Time.now
hand1 = cards[0..12]
hand2 = cards[13..25]
hand3 = cards[26..38]
hand4 = cards[39..51]
p_sort = 0
[hand1, hand2, hand3, hand4].each { |hand|
p_sort += 1 if perfect_hand_by_sort?(hand)
}
dt_sort = Time.now - t0
t_by_sort += dt_sort
p_scan = 0
t0 = Time.now
[1,2,3,4].each { |handno|
p_scan += 1 if perfect_hand_by_scan?(cards, handno)
}
dt_scan = Time.now - t0
t_by_scan += dt_scan
p_minmax = 0
t0 = Time.now
[1,2,3,4].each { |handno|
p_minmax += 1 if perfect_hand_by_minmax?(cards, handno)
}
dt_minmax = Time.now - t0
t_by_minmax += dt_minmax
perfect_hands += p_sort
end
puts "After #{TRIALS} trials, we got #{perfect_hands}"
puts "Times: "
puts " sort: #{t_by_sort}"
puts " scan: #{t_by_scan}"
puts " minmax: #{t_by_minmax}"
Problem 9 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 9: What are the odds of winning in Craps (using only Pass Line bets)?
#!/usr/bin/env ruby
# Modelling craps is weird.
# We roll two dice and sum them.
#
# if they're 7 or 11 => immediate win
# 2,3,12 => lose
# else, roll is the "point"
# player keeps rolling until:
# point => win
# 7 => lose
#
# So, what's the chance to win?
TRIALS=500000
def roll
d1 = 1 + rand(6)
d2 = 1 + rand(6)
return d1+d2
end
def win()
first_roll = roll()
return true if (first_roll == 7 || first_roll == 11)
return false if ([2,3,12].include?(first_roll))
point = first_roll
while (true) do
next_roll = roll()
return false if 7 == next_roll
return true if next_roll == point
end
end
wins = 0
TRIALS.times do
wins += 1 if win()
end
puts "After #{TRIALS} games, won #{wins}"
puts "P(win) = #{wins.to_f/TRIALS.to_f}"
Problem 10 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 10: An urn drawing game: an urn contains 10 black and 10 white balls. You call "black" or "white," and a ball is drawn. If it matches your call, you win $10. a) What's the most you will pay to play? b) Change the game such that the black/white mix is unknown to you. What would you pay if you were only allowed to play the game one time?
#!/usr/bin/env ruby
TRIALS=10000
# For the fair game, we can just try it a bunch of times.
# The game is symmetric about color choice, so we can
# ignore the details and just use the 50/50 odds.
wins = 0
TRIALS.times do
wins +=1 if (rand() < 0.5)
end
payout = 10.0 * wins
puts "Out of #{TRIALS} fair plays, we won #{wins}, $#{payout}."
puts "Therefore, the maximum amount to pay would be "
puts " $#{payout/TRIALS.to_f}"
# Now, for the weird game. Your friend can vary the
# ratio of balls in the urn, so we need to try a lot
#of games at different probabilities. We can do this
# stratified or randomly. Let's do randomly first.
wins = 0
TRIALS.times do
wins +=1 if (rand() < rand())
end
payout = 10.0 * wins
puts "Out of #{TRIALS} randomized friend-chosen plays"
puts " we won #{wins}, $#{payout}. Therefore, the "
puts " maximum amount to pay would be "
puts "$#{payout/TRIALS.to_f}"
# And now stratified, to see if it makes any difference
wins = 0
TRIALS.times do |i|
wins +=1 if (rand() < i.to_f/TRIALS.to_f)
end
payout = 10.0 * wins
puts "Out of #{TRIALS} stratified friend-chosen plays,"
puts " we won #{wins}, $#{payout}. Therefore, the "
puts " maximum amount to pay would be "
puts " $#{payout/TRIALS.to_f}"
Problem 15 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 15: 8 boys and 7 girls are seated at random in a row of 15 seats. On the average, how many pairs of adjacent seats are taken by cootie-transmitting couples? (boys next to girls)
#!/usr/bin/env ruby
TRIALS=10000
# The original problem is heteronormative bachelors
# and models. I don't like that for various reasons,
# and cooties is funnier anyway.
#
# But, my code represents the "bachelors and models in
# marriagable couples" model of the original problem.
#
# We need a pool of bachelors and models:
N_BACHELORS=8
N_MODELS=7
class Pool
def initialize
@b = N_BACHELORS
@m = N_MODELS
end
def draw
left = @b+@m
return nil if left == 0
r = rand(left)
if r > @b
@m -= 1
return "M"
else
@b -= 1
return "B"
end
end
end
n_couples = 0
TRIALS.times do
p = Pool.new
last_d = nil
while d = p.draw()
n_couples += 1 if last_d && last_d != d
last_d = d
end
end
puts "After #{TRIALS}, got #{n_couples}."
puts "So, #{n_couples/TRIALS.to_f} on average"
Problem 16 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 16: A tennis tournament of 8 players has the initial ladder chosen at random. Assuming the best player beats everyone, and the second-best player beats everyone else, what's the chance that the second-best player wins runner-up?
#!/usr/bin/env ruby
TRIALS=100000
# Here we get into some slightly gnarly modelling.
#
# We'll model the current stage of competition as an array
# of ranks:
PLAYERS=[1,2, 3,3, 3,3, 3,3]
# Then, we run a round of the game as follows:
def round(seedings)
next_round = []
i = 0
while i < seedings.length
a = seedings[i]
b = seedings[i+1]
winner = a < b ? a : b # ranks, so lower wins
next_round.push(winner)
i += 2
end
return next_round
end
def second_runner_up(first_round)
quarters = round(first_round)
semis = round(quarters)
return semis.include?(2)
end
n_second_runner_up = 0
TRIALS.times do
srand # There's a fun story here... guess what it is.
seeding = PLAYERS.shuffle
n_second_runner_up += 1 if second_runner_up(seeding)
end
puts "Out of #{TRIALS} brackets, #2 was runner up "
puts " #{n_second_runner_up} times. "
puts " P = #{n_second_runner_up/TRIALS.to_f}"
Problem 17 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 17: In a jousting tournament, there are two twins: Balan and Balin. If there are 8 knights total, all equally matched, and the initial ladder is randomly assigned, what's the probability of the twins jousting against each other?
#!/usr/bin/env ruby
TRIALS=100000
# Array of Knighs: BalIn and BalAn, and some schmucks.
KNIGHTS=["I","A","X","X","X","X","X","X"]
# The knights are equally matched, so we need to randomize
# the selection of who wins each round. Otherwise, this is
# very similar to problem16.rb.
# We're going to abuse the weak typing of ruby here. Don't
# do this for real.
#
# Returns an array if we just continue on, returns nil if the
# brothers meet.
#
def round(seedings)
next_round = []
# Walk the seedings pairwise, looking for "A" and "I"
# Note that order doesn't matter.
i = 0
while (i < seedings.length())
return nil if (["A","I"] & seedings[i..i+1]).length == 2 # rubylicious.
next_round.push( seedings[i + rand(2)] )
i += 2
end
return next_round
end
def tourney()
first_round = KNIGHTS.shuffle
quarters = round(first_round)
return true if !quarters
semi = round(quarters)
return true if !semi
final = round(semi)
return true if !final
return false # tourney over
end
srand() # Our first call to shuffle is before our
# first call to rand()...
n_match=0
TRIALS.times do
n_match += 1 if tourney()
end
puts "In #{TRIALS} trials, the brothers met #{n_match} times."
puts "P=#{n_match/TRIALS.to_f}"
Problem 18 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 18: What's the probability of getting 50 tails in 100 toin cosses?
#!/usr/bin/env ruby
# This is pretty straightforward:
# we just count the # of 50s we get.
TRIALS=1000000
n_fifties = 0
TRIALS.times do
heads = 0
100.times { heads += rand(2) }
n_fifties += 1 if 50 == heads
end
puts "Out of #{TRIALS} trials"
puts " we got #{n_fifties} even splits, so "
puts " P=#{n_fifties/TRIALS.to_f}"
Problem 19 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 19: Which is most likely: a) at least 1 six when rolling 6 dice, b) at least 2 sixes in 12 dice rolled, or c) at least 3 sixes of 18 dice?
#!/usr/bin/env ruby
TRIALS=10000
n_ones = 0
n_twos = 0
n_threes = 0
def roll
return 1+rand(6)
end
TRIALS.times do
six = (1..6).to_a.map { roll() }
twelve = (1..12).to_a.map { roll() }
eighteen = (1..18).to_a.map { roll() }
n_ones += 1 if six.select { |r| r == 6 }.length >= 1
n_twos += 1 if twelve.select { |r| r == 6 }.length >= 2
n_threes += 1 if eighteen.select { |r| r == 6 }.length >= 3
end
puts "After #{TRIALS}, got #{n_ones} / #{n_twos} / #{n_threes}"
Problem 20 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 20: in a stepwise, three cornered duel between Adam, Bryan, and Costello, what should Adam's strategy be? Adam shoots first, then Bryan, then Costello, then Adam, and so on, until only one stands. The probability of Adam hitting his target is 0.3, Bryan is a perfect shot, and C is 50/50.
Note that the code below doesn't actually answer the question the same way as Mosteller. He admits a possibility our code below doesn't, which, personally, I feel is a breach of honor in such an esteemed tradition as settling arguments by the well-reasoned method of shooting at each other.
#!/usr/bin/env ruby
# I don't like the answer in the book here;
# code of honor is rough.
class Duel
P_A = 0.3
P_B = 1.0
P_C = 0.5
def initialize()
@a_live = true
@b_live = true
@c_live = true
@a_moves = []
end
attr_reader :a_moves
def over?()
!@a_live || !(@b_live || @c_live)
end
def a_live?()
return @a_live
end
def run_round(moves)
@a_moves.push(moves[0])
if rand() < P_A
if moves[0] == 0 && @b_live
@b_live = false
else
@c_live = false
end
end
if @b_live && rand() < P_B
if moves[1] == 0 || !@c_live
@a_live = false
else
@c_live = false
end
end
if @c_live && rand() < P_C
if moves[1] == 0 && @b_live
@b_live = false
else
@a_live = false
end
end
end
end
move_space = []
8.times { |move_mask|
move_a = move_mask & 1 == 0 ? 1 : 0
move_b = move_mask & 2 == 0 ? 1 : 0
move_c = move_mask & 4 == 0 ? 1 : 0
move_space.push( [move_a, move_b, move_c] )
}
a_lives = Hash.new { |h,k| h[k] = 0 }
a_attempts = Hash.new { |h,k| h[k] = 0 }
two_move_space = []
al1 = [0,0]
aa1 = [0,0]
move_space.each { |m1|
a1 = m1[0]
move_space.each { |m2|
100.times {
d = Duel.new
d.run_round(m1)
d.run_round(m2)
while (!d.over?)
d.run_round(m2) # doesn't matter, we only have c left if we're still in the game
end
aa1[a1] += 1
al1[a1] += 1 if d.a_live?
a_lives[ d.a_moves ] += 1 if d.a_live?
a_attempts[ d.a_moves ] += 1
}
}
}
puts a_lives.inspect
puts a_attempts.inspect
lives = [0,0]
attempts = [0,0]
a_lives.each { |k,v|
first_move = k[0]
att = a_attempts[k]
lives[first_move] += v
attempts[first_move] += att
}
puts lives.inspect
puts attempts.inspect
puts
puts aa1.inspect
puts al1.inspect
puts [ al1[0]/aa1[0].to_f, al1[1]/aa1[1].to_f ].inspect
Problem 35 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 35: After imbibing a few too many, a drunk man finds himself stumbling at the edge of a cliff. One step to the left will send him over, so he's trying to move right as well as he can. Except that he's drunk, so he only moves to the right with probability 2/3. What are his chances of living into sobriety?
#!/usr/bin/env ruby
TRIALS=10000
POS_START=1 # Initial position
# Probability of a step in the positive direction
P_POSITIVE=0.666667
# What position he should be at to "get away"
#
# Note that this doesn't actually guarantee that our friend
# lives to the morrow, but it's close enough.
#
GOT_AWAY = 100000
MAX_STEPS = 1000000 # How long do we wait?
wins=losses=0
TRIALS.times do
pos = POS_START
steps = 0
while (pos > 0 && pos < GOT_AWAY && steps < MAX_STEPS)
pos += rand() < P_POSITIVE ? 1 : -1
steps += 1
end
losses += 1 if (pos == 0)
wins += 1 if (pos > 0)
end
print "Wins: #{wins}, losses: #{losses}"
Problem 37 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 37: We need $40 to get on the bus home from Vegas tomorrow, but are down to $20. The plan is to play evens in roulette (2:1 payout, 18/38 probability of winning), but we're missing one detail. Do we bet it all one time and walk away with $0 or $40, or do we bet it a dollar at a time?
#!/usr/bin/env ruby
# Bold play or cautious play? We need $40 to get on the
# bus home from Vegas tomorrow, but are down to $20. Do
# we bet it all at once on evens in roulette, or do we
# bet it a dollar at a time?
TRIALS=10000
P_WIN = 18.0/38.0 # 18 evens on a 38-slot roulette wheel
def play_bold()
return rand() < P_WIN
end
def play_cautious()
bank = 20
plays_remaining = 10000 # limit how long we can play
while (bank < 40 && bank > 0 && plays_remaining > 0)
if (rand() < P_WIN)
bank += 1
else
bank -= 1
end
plays_remaining -= 1
end
return bank == 40
end
win_bold = 0
win_cautious = 0
TRIALS.times {
win_bold += 1 if play_bold()
win_cautious += 1 if play_cautious()
}
puts "After #{TRIALS} times, wins:"
puts " bold: #{win_bold}"
puts " cautious: #{win_cautious}"
Problem 39 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 39: a set of glass rods, with one end blue-dotted and the other red-dotted, are dropped. Many break into 3 parts. What is the average length of the blue-dotted pieces of these broken-in-three rods?
#!/usr/bin/env ruby
TRIALS=100000
# Much like 42 and 43, this is pretty trivial
cum_blue = 0.0
TRIALS.times {
# Choose two breaking points
b1 = rand()
b2 = rand()
# put them on a stick reasonably...
if (b1 > b2)
x = b1
b1 = b2
b2 = x
end
# and let's say the blue is always on the left
# that is, the length of the blue segment is b1
cum_blue += b1
}
puts "After #{TRIALS} trials, average length:"
puts cum_blue/TRIALS
Problem 42 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 42: If a stick is broken in two at random, what is the average length of the shorter half? What is the average ratio of the shorter side to the longer side?
#!/usr/bin/env ruby
TRIALS = 1000000
# This is pretty trivial: we just choose a breakpoint somewhere
# along a unit-length stick and keep track of the total
# length of the short ends, and the ratio for (b)
cum_short_len = 0.0
cum_ratio = 0.0
TRIALS.times {
breakpoint = rand()
breakpoint = 1.0 - breakpoint if breakpoint > 0.5
cum_short_len += breakpoint
cum_ratio += breakpoint / (1.0 - breakpoint)
}
mean = cum_short_len / TRIALS
mean_ratio = cum_ratio / TRIALS
puts "After #{TRIALS} attempts, mean(short_len) = #{mean}"
puts " mean(ratio) = #{mean_ratio}"
Problem 43 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 43: A bar is broken at random in two places. Find the average length of the shortest, middle-est, and longest pieces.
#!/usr/bin/env ruby
TRIALS=100000
# Much like 42, this is pretty trivial stuff.
cum_short = 0.0
cum_med = 0.0
cum_long = 0.0
TRIALS.times {
b1 = rand()
b2 = rand()
# put them on a stick reasonably...
if (b1 > b2)
x = b1
b1 = b2
b2 = x
end
ls = [b1, b2-b1, 1.0-b2].sort
cum_short += ls[0]
cum_med += ls[1]
cum_long += ls[2]
}
m_short = cum_short/TRIALS
m_med = cum_med/TRIALS
m_long = cum_long/TRIALS
puts "After #{TRIALS} trials, average lengths:"
puts " #{m_short}, #{m_med}, #{m_long}"
Problem 45 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 45: Let us shuffle two decks of cards and lay them out in two lines, one over the other. On average, how many cards will "line up" with themselves (ie: the 3 of clubs over the 3 of clubs)?
#!/usr/bin/env ruby
TRIALS=10000
deck = (0..51).to_a
# We shuffle the deck, then shuffle another. We then walk the
# two arrays together, looking for matches.
match = 0
TRIALS.times do
upper = deck.shuffle
lower = deck.shuffle
upper.each_with_index do |u_i, i|
l_i = lower[i]
match += 1 if l_i == u_i
end
end
puts "Out of #{TRIALS}, #{match} cards matched."
puts "That gives us an expected number of matches as: "
puts " #{match.to_f/TRIALS.to_f}"
Problem 46 of Monte Carlo solutions to Fifty Challenging Problems...
Problem 46: As in 45, but this time: what is the probability of r matches?
#!/usr/bin/env ruby
TRIALS=1000000
deck = (0..51).to_a
# We shuffle the deck, then shuffle another.
# We then walk the two arrays side-by-side, checking
# for matches
# a count of how many times we got a given number of matches
matches = deck.map { 0 }
TRIALS.times do
upper = deck.shuffle
lower = deck.shuffle
match = 0
upper.each_with_index do |u_i, i|
l_i = lower[i]
match += 1 if l_i == u_i
end
matches[match] += 1
end
puts "Out of #{TRIALS}, we got the following distribution:"
matches.each_with_index { |m_i, i| puts "\t#{i}\t#{m_i}" }
This was brought to you by Josh Myer. He has other fun things at his homepage.