## 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.