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.