5.3. Monte Carlo

using PyPlot
using Random

5.3.1. Monte Carlo simulations

A common technique to estimate probabilities and statistical quantities is Monte Carlo simulation. The computer runs many simulations that depend on random numbers and records the outcomes. The desired quantity can then be estimated by simple fractions.

As an example, the code below simulates rolling a die repeatedly until it rolls a 6, and recording the number of rolls that it took. The average number of rolls required can then be estimated.

ntrials = 10000            # Number of experiments
x = Int64[]                # Array of outcomes
sumx = 0
for i = 1:ntrials
    nrolls = 0
    found_a_six = false
    while !found_a_six
        nrolls += 1
        roll = rand(1:6)
        if roll == 6
            found_a_six = true
        end
    end
    push!(x, nrolls)
    sumx += nrolls
end

plt.hist(x, 50);
xlabel("Number of rolls")
ylabel("Count")
average = sumx / ntrials

5.9805

5.3.1.1. Card games, permutations

To simulate card games, we need a way to represent the cards. There are 13 ranks and 4 suits, so we could simply use an integer between 1 and 52 and let 1 to 13 represent Ace, 2, 3, …, King of clubs, 14 to 26 all the diamonds, etc. We can use integer division and remainder to find the rank and suit:

function card_rank(card)
    (card - 1) % 13 + 1
end

function card_suit(card)
    (card - 1) ÷ 13 + 1
end

cards = 1:52      # All cards
[card_rank.(cards) card_suit.(cards)]

52×2 Matrix{Int64}:
  1  1
  2  1
  3  1
  4  1
  5  1
  6  1
  7  1
  8  1
  9  1
 10  1
 11  1
 12  1
 13  1
  ⋮  
  2  4
  3  4
  4  4
  5  4
  6  4
  7  4
  8  4
  9  4
 10  4
 11  4
 12  4
 13  4

We can now draw a random card by generating a random integer between 1 and 52. However, to deal e.g. a poker hand of 5 cards, we need to make sure we do not draw the same card twice. One way to do this is with Julia’s randperm(n) function, which generates a random permutation of the integers 1 up to n. We can use this to represent a shuffled deck of cards, and use the first 5 cards to deal a poker hand:

cards = randperm(52)
hand = cards[1:5]
[card_rank.(hand) card_suit.(hand)]

5×2 Matrix{Int64}:
  4  1
  1  3
 11  1
  3  4
 10  3

5.3.1.2. Example: Probability of flush

A flush means that all cards have the same suit. To estimate the probabilty that a random poker hand is a flush, we run Monte Carlo simulation:

ntrials = 100000
nflush = 0
for itrial = 1:ntrials
    cards = randperm(52)
    hand = cards[1:5]

    suits = card_suit.(hand)
    same_suit = true
    for i = 2:5
        if suits[i] ≠ suits[1]
            same_suit = false
            break
        end
    end
    if same_suit
        nflush += 1
    end
end

approx_probability = nflush / ntrials

0.00205

This problem can also be solved using combinatorial techniques, which we can use to check the approximation:

exact_probability = 4 * binomial(13,5) / binomial(52,5)

0.0019807923169267707

5.3.2. Estimate \(\pi\) by throwing darts

Random number can be used to approximate areas, or more generally, to estimate integrals. To illustrate this, consider a 2-by-2 square with a unit disk inside. Imagine throwing \(n\) darts on the square with equal probability everywhere. Then the fraction of darts that hit inside the circle should approximate the ratio between the areas, that is:

\[ \frac{\mathrm{hits}}{n} \approx \frac{\pi r^2}{(2r)^2} \]

where the radius \(r=1\). This means we can estimate \(\pi\) as

\[ \pi \approx 4 \frac{\mathrm{hits}}{n} \]

n = 10000               # Number of darts
x = 2rand(n) .- 1       # Coordinates
y = 2rand(n) .- 1

plot([-1,1,1,-1,-1], [-1,-1,1,1,-1])  # Draw square
theta = 2π*(0:100)./100               # Draw circle
plot(cos.(theta), sin.(theta))   
plot(x, y, linestyle="None", marker=".", markersize=1)   # Plot dart points
axis("equal")

# Determine if points are inside the circle (a "hit")
hits = 0
for i = 1:n
    if x[i]^2 + y[i]^2 ≤ 1
        hits += 1
    end
end

approx_pi = 4hits / n

3.1424

5.3.3. Example: Random walk

A random walk can be described by the following algorithm:

• Consider a 2-d array of square cells

• Start at the center cell

• At each turn, randomly choose a direction (up/down/left/right) and move to that neighboring cell

• Continue until reaching an edge

function random_walk(n)
    x = [0]
    y = [0]
    while abs(x[end]) < n && abs(y[end]) < n
        if rand() < 0.5
            if rand() < 0.5  # Up
                push!(x, x[end])
                push!(y, y[end] + 1)
            else             # Down
                push!(x, x[end])
                push!(y, y[end] - 1)
            end
        else
            if rand() < 0.5  # Right
                push!(x, x[end] + 1)
                push!(y, y[end])
            else             # Left    
                push!(x, x[end] - 1)
                push!(y, y[end])
            end
        end
    end
    x,y
end    

random_walk (generic function with 1 method)

n = 10
x,y = random_walk(n)

plot(x, y, marker=".", markersize=8)     # Draw dots at each point in x,y
grid(true)
axis([-n,n,-n,n]);