Archive for the ‘Uncategorized’ Category

Monte Carlo: part 23

Friday, September 3rd, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 22

Tuesday, August 31st, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 21

Saturday, August 28th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 20

Wednesday, August 25th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 19

Sunday, August 22nd, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}" }

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 18

Thursday, August 19th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 17

Monday, August 16th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 16

Friday, August 13th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 15

Tuesday, August 10th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.

Monte Carlo: part 14

Saturday, August 7th, 2010

(This is another part of the Fifty Problems series, a set of example applications of Monte Carlo methods. In each post, I present source code which answers a probabilistic question using simulated models of the underlying system.)

Problem 18: What’s the probability of getting 50 tails in 100 toin cosses?

#!/usr/bin/env ruby

# This is pretty straightforward modeling:
#   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}"

I’ve been coding my way through Fifty Challenging Problems in Probabilities with Solutions. This post is a part of the Fifty Challenging Problems series. Read all of the Fifty Problems posts by clicking here.