Introduction

Monte Carlo methods are named after the city in Monaco which is known for it’s casinos. Credit for inventing the method is often given to Stnislaw Ulam, a mathematician who worked closely with Joh von Neumann on the Manhattan Project for the United States during WWII.

The basic idea of Monte Carlo is to approximate probabilities or expected values, by simulating a large number of samples from the distribution of interest. Assume that we have a function \(f(x)\) which represents either a PMF (discrete) or PDF (continuous). If we can obtain \(M\) indepdent and identically distributed draws, \(x_1, x_2, \cdots x_M\) from the distribution, then we can approximate \(f\) by \[\frac{1}{M}\sum_{i=1}^M 1(x_i=x) \stackrel{}{\rightarrow} f(x)\]

Note that \(1()\) is the indicator function and is defined as \[1(x_i=x) = \begin{cases} 1, & x_i = x \\ 0, & x_i \neq x \end{cases}\]

Let’s not get too bogged down in the technical details, Monte Carlo is actually quite simple.

Example 1 - Estimating \(\pi\) with Monte Carlo

Consider a square centered at the origin having sides of length \(2\), and a circle with radius \(1\) placed inside this square.

We can choose a point from this square uniformly at random by doing \[X \sim U(-1, 1)\] \[Y \sim U(-1, 1)\] with \(X\) and \(Y\) independent. What is the probability that this point lands inside the circle?

This problem is actually not so difficult analytically, but we can use a Monte Carlo approach to approximate an answer.

  1. We draw a large number \((M)\) of points from the rectangle.
  2. Count how many of the points are inside the circle and divide by \(M\).
#Choose M
M <- 100
#Draw M points from rectangle
x <- runif(M, -1, 1)
y <- runif(M, -1, 1)
#Calculate radius of each point
r <- sqrt(x^2 + y^2)
#Determine which points are in circle
in_circle <- (r < 1)

Since there are 78 points inside the circle, we can estimate the probability with

sum(in_circle)/M
## [1] 0.78

Of course, this probability approximation should get better if we increase \(M\). With \(M=10,000\) we get

#Repeat for M=10,000
sum(in_circle)/M
## [1] 0.7857

It is important to note, for this problem we didn’t really need Monte Carlo. In fact with uniform sampling, the probability is just \[P(\text{in circle}) = \frac{\text{area of circle}}{\text{area of square}} = \frac{\pi\cdot 1^2}{2^2} = \frac{\pi}{4} \approx 0.785\]

Rearranging the terms gives \[\pi = 4\cdot P(\text{in circle})\]

So if we plug in our approximation for \(P(\text{in circle})\) to this equation, we obtain a Monte Carlo estimate for \(\pi\).

4*sum(in_circle)/M
## [1] 3.1428

As a final illustration, lets see how the estimate looks as a function of \(M\).

Example 2 - Monte Carlo for Two Variables

Consider Newton’s second law of motion \[F= ma\] which states that Force = mass \(\times\) acceleration.

Suppose a random force \(F \sim N(100, 4^2)\) (Newtons) is applied to an object with random mass \(M \sim Gamma(10,1)\) (Kg). What is the mean and variance of the resulting random acceleration \(A\)?

This problem can be very easily approximated with Monte Carlo.

#Choose M
M <- 10000
#Simulate Force and Mass
Force <- rnorm(M, mean=100, sd=4)
Mass  <- rgamma(M, shape=10, rate=1)
#Calculate random Acceleration
A <- Force/Mass
#Create histogram of A
hist(A, breaks=20, main='Monte Carlo Distribution of Random Acceleration')

Now we can approximate the mean and variance with the sample mean and sample variance.

mean(A)
## [1] 11.15101
var(A)
## [1] 14.92043

Example 3 - A Monte Carlo Solution to the Labyrinth Problem

Consider the challenge problem from Homework 3. We can re-do this problem with Monte Carlo (and actually obtain a better answer for part b). Simulating this problem is somewhat trickier… I’ll try my best to comment what I’m doing.

#Number of tunnels 
n <- 7
#Initally the time is 0
X <- 0
#Set the still trapped flag to TRUE
still_trapped <- TRUE

while(still_trapped){
  #Choose a tunnel at random
  tunnel <- sample(n, 1)
  #Add tunnel time to X
  X <- X + tunnel
  #Check to see if we escaped
  if(tunnel == n){
    still_trapped <- F
  }
}

For this particular simulation, it takes us X= 22 hours to escape. For simplicity going forward, I’m going to define this as a function in R. It will take a single argument, n the number of tunnels, and it will return the time it takes to escape for a single simulation.

time2escape <- function(n){
  X <- 0
  still_trapped <- TRUE
  while(still_trapped){
    tunnel <- sample(n, 1)
    X <- X + tunnel
    if(tunnel == n){
      still_trapped <- F
    }
  }
  return(X)
}

Now we can estimate the distribution of \(X\) using Monte Carlo! We start by fixing \(n = 7\).

#Choose M
M <- 10000
#Make X an initially empty vector of length M
X <- rep(NA, M)
#Do the simulation M times
for(i in 1:M){
  X[i] <- time2escape(n=7)
}
hist(X, main='Monte Carlo Distribution for Escape Time (n=7)', breaks=20)

So to answer the question How many hours will it take you to escape on average?, we simply check

mean(X)
## [1] 28.2582

We can repeat this process for different tunnel lengths and plot the answer as a function of \(n\). This plot demonstrates that, on average, we will not escape for more than \(n=6\) tunnels.

For part b) in the homework problem, finding specific probabilities was a challenging (impossible?) task analytically. The best we were able to do was find a lower bound using Markov’s Inequality. As long as we are okay with getting an approximate answer, we can use Monte Carlo.

We can do this with very minor additions to the old code.

#Choose M
M <- 10000
#Make X an initially empty vector of length M
X <- escape <- rep(NA, M)
#Do the simulation M times
for(i in 1:M){
  X[i] <- time2escape(n=7)
  #Did we escape?
  escape[i] <- (X[i] <= 24) #Returns either TRUE or FALSE
}
hist(X, main='Monte Carlo Distribution for Escape Time (n=7)', breaks=20)
abline(v=24, lty=3, col='red')

We can estimate the probability with

sum(escape)/M
## [1] 0.571

We have learned something interesting from this analysis. Note that for \(n=7\) tunnels, we said that we will fail to escape on average. However, the probability of escaping is 0.571, which is better than \(1/2\).

Finally, as before, we can plot this is a function of \(n\).

Assignment

Problem One

Let \(Z \sim N(0,1)\). Use a Monte Carlo Algorithm to approximate

  • \(P(|Z| < 1)\)
  • \(P(|Z| < 2)\)
  • \(P(|Z| < 3)\)

Note: This is known as the \(68-95-99.7\) rule.

Problem Two

The Logistic Equation is sometimes used as model of population growth. Suppose that on a small island, there are \(X_0\) bunnies at time \(t=0\). The island has a carrying capacity of \(K\) bunnies (the maximum number of bunnies that can survive on the island). The population has a growth rate (birth rate - death rate) of \(R\). At time \(t\), the Logistic Model says that the number of rabbits on the island will be

\[X_t = \frac{ \exp(Rt)}{1/X_0 + (\exp(Rt)-1)/K}\] Professor Halfbrain is studying these rabbits. Although he has tried to determine \(X_0\), \(K\) and \(R\), he is uncertain. Therefore he supplements his model with \[X_0 \sim Poisson(30)\] \[K' \sim N(200, 5^2) \quad and \quad K = \lfloor K' \rfloor\] \[R \sim Beta(2, 20)\]

Use a Monte Carlo algorithm to approximate the expected value and variance of \(X_t\) after \(t=12\) months. Also approximate \(P(X_t > 75)\).

Helpful Functions

  • The rnorm(n, mu, sigma) function generates n random variables from the Normal distribution with mean mu and standard deviation sigma.
  • The abs(x) function is the "absolute value function". Examples: abs(3.14) = 3.14 and abs(-3.14) = 3.14.
  • The rpois(n, lambda) function generates n random variables from the Poisson distribution with mean lambda.
  • The rbeta(n, alpha, beta) function generates n random variables from the Beta distribution with parameters alpha and beta
  • The floor(x) function is the "floor function". Examples: floor(3.5) = 3, floor(3) = 3, floor(12.001) = 3, floor(12.99) = 12