Monte Carlo Overview
Simulation in R
Simulation of Roulette
Consider the casino game Roulette.
We can use simulation to evaluate gambling strategies.
Roulette Simulation in R
RouletteSpin <- function(num.spins){
# function to simulate roulette spins
# ARGS: number of spins
# RETURNS: result of each spin
outcomes <- data.frame(number = c('00','0',
as.character(1:36)),
color=c('green','green','red','black','red','black',
'red','black','red','black','red','black',
'black','red','black','red','black',
'red','black','red','red','black',
'red','black','red','black','red',
'black','red','black','black','red',
'black','red','black','red','black','red'))
return(outcomes[sample(38,num.spins,replace=T),])
}
kable(RouletteSpin(2), row.names=F)| number | color |
|---|---|
| 34 | red |
| 13 | black |
Monte Carlo Procedures
Motivation for Monte Carlo procedures
Some probabilities can easily be calculated either intuitively or using pen and paper; however, as we have seen we can also simulate procedures to get an approximate answer.
Consider poker, where players are dealt a hand of five cards from a deck of 52 cards. What is the probability of being dealt a full house?
Full House Probability Calculation
Could we analytically compute this probability? Yes Is it an easy calculation? not necessarily. So consider a (Monte Carlo) simulation.
Full House Probability Calculation
Could we analytically compute this probability? Yes Is it an easy calculation, not necessarily. So consider a (Monte Carlo) simulation.
Full House Probability Calculation
IsFullHouse <- function(hand){
#determines whether a hand of 5 cards is a full house
#ARGS: data frame of 5 cards
#RETURNS: TRUE or FALSE
cards <- hand[,2]
num.unique <- length(unique(cards))
num.appearances <- aggregate(rep(1,5),
list(cards),sum)
max.appearance <- max(num.appearances[,2])
ifelse(num.unique == 2 & max.appearance ==3,
return(TRUE),return(FALSE))
}Full House Probability Calculation
num.sims <- 1e5
sim.hands <- replicate(num.sims,DealPoker()
,simplify=F)
results <- lapply(sim.hands,IsFullHouse)
prob.full.house <- sum(unlist(results))/num.simsAnalytically the probability of getting a full house can be calculated as approximately 0.00144, with our simulation procedure we get 0.00171.
Closing Thoughts on Monte Carlo
A few facts about Monte Carlo procedures:
- They return a random result due to randomness in the sampling procedures.
- The run-time (or number of iterations) is fixed and typically specified.
- Mathematically, Monte Carlo procedures are computing an integral or enumerating a sum.
- They take advantage of the law of large numbers
- They were given the code name Monte Carlo in reference to the Monte Carlo Casino in Monaco.
Exercise: Probability of Red, Green, and Black
Calculate/Derive the probability of landing on green, red, and black.
How can the
RouletteSpin()function be used to compute or approximate these probabilities?
Solution: Probability of Red, Green, and Black
In this situation, it is easy to compute the probabilities of each color analytically. However, consider simulating this process many times to estimate these probabilities.
num.sims <- 1000
spins <- RouletteSpin(num.sims)
p.red <- sum(spins[,2] == 'red') / num.sims
p.black <- sum(spins[,2] == 'black') / num.sims
p.green <- sum(spins[,2] == 'green') / num.simsAnalytically \(P[red] = \frac{18}{38} =\) 0.4737, this is estimated as 0.477. Similarly, \(P[black] = \frac{18}{38} =\) 0.4737, this is estimated as 0.472 and \(P[green] = \frac{2}{38} =\) 0.0526, this is estimated as 0.051
Exercise: Simulation Questions - Part 2
Now what happens if we:
run the simulation again with the same number of trials?
run the simulation with more trials, say 1 million?
Solution: Simulation Questions - Part 2
Now what happens if we:
- run the simulation again with the same number of trials?
num.sims <- 1000
spins <- RouletteSpin(num.sims)
p.red <- sum(spins[,2] == 'red') / num.sims
p.black <- sum(spins[,2] == 'black') / num.sims
p.green <- sum(spins[,2] == 'green') / num.simsThe simulated results are different Analytically \(P[red] = \frac{18}{38} =\) 0.4737, this is estimated as 0.456. Similarly, \(P[black] = \frac{18}{38} =\) 0.4737, this is estimated as 0.488 and \(P[green] = \frac{2}{38} =\) 0.0526, this is estimated as 0.056
Solution: Simulation Questions - Part 2
- run the simulation with more trials, say 1 million?
num.sims <- 1000000
spins <- RouletteSpin(num.sims)
p.red <- sum(spins[,2] == 'red') / num.sims
p.black <- sum(spins[,2] == 'black') / num.sims
p.green <- sum(spins[,2] == 'green') / num.simsAnalytically \(P[red] = \frac{18}{38} =\) 0.4737, this is estimated as 0.4742. Similarly, \(P[black] = \frac{18}{38} =\) 0.4737, this is estimated as 0.4731 and \(P[green] = \frac{2}{38} =\) 0.0526, this is estimated as 0.0528
Yahtzee
Implement a Monte Carlo procedure to estimate the probability of rolling a Yahtzee (5 dice all with the same value) with a single roll of the dice.
Yahtzee: Part II
Suppose that you have currently rolled 1, 3, 3, 6, 4. If you keep the pair of 3’s and re-roll the other three dice, what is the probability of rolling a Yahtzee in the next roll?
(For those of you that play Yahtzee, you know that you get a total of three rolls. You could also simulated a third roll too.)
Black-Jack
Implement a procedure to estimate the probability of being dealt black-jack (2 cards that total 21).
Black-Jack: Part II
Assume you have been dealt the 10 of hearts and the 4 of diamonds and the dealer has a 6 showing. By standard rules, you win if your card total is higher than the dealer, but less than or equal to 21. You have the option to hit which means adding more cards to your total, but if you go over 21 you bust and automatically lose.
Describe in words (or pseudeocode) how you would simulate this process to estimate the probability of winning by hitting and staying.