---
title: "Monte Carlo Overview"
output:
  revealjs::revealjs_presentation:
    theme: white
    center: true
    transition: none
    incremental: true
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache=TRUE)
library(knitr)
library(tidyverse)
set.seed(09222021)
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),])
}
```


# Simulation in R

## Simulation of Roulette
Consider the casino game Roulette.

![Roulette Wheel](roulette.jpg)

We can use simulation to evaluate gambling strategies.

## Roulette Simulation in R
```{r roulette}
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)
```


# 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.
```{r}
DealPoker <- function(){
  # Function returns a hand of 5 cards
  
  
  
  
}
  
```

## Full House Probability Calculation
Could we analytically compute this probability? **Yes** Is it an easy calculation, not necessarily. So consider a (Monte Carlo) simulation.
\footnotesize
```{r}
DealPoker <- function(){
  # Function returns a hand of 5 cards
  deck <- data.frame( suit = rep(c("H", "S","C","D"),each=13),
          card = rep(c('A',2:10,'J',"Q",'K'),4) )  
  return(deck[sample(52,5),])
}
```

## Full House Probability Calculation
\footnotesize
```{r}
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
```{r}
num.sims <- 1e5
sim.hands <- replicate(num.sims,DealPoker()
                       ,simplify=F)
results <- lapply(sim.hands,IsFullHouse)
prob.full.house <- sum(unlist(results))/num.sims
```
Analytically the probability of getting a full house *can* be calculated as approximately 0.00144, with our simulation procedure we get `r prob.full.house`.

## 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
1. Calculate/Derive the probability of landing on green, red, and black.

2. 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.
```{r rouletteSpins}
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.sims
```
Analytically $P[red] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.red,4)`. Similarly, $P[black] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.black,4)` and $P[green] = \frac{2}{38} =$ `r round(2/38,4)`, this is estimated as `r round(p.green,4)`

## Exercise: Simulation Questions - Part 2 

Now what happens if we:

1. run the simulation again with the same number of trials?

2. run the simulation with more trials, say 1 million?

## Solution: Simulation Questions - Part 2

Now what happens if we:

1. run the simulation again with the same number of trials?
\footnotesize
```{r }
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.sims
```
The simulated results are different
Analytically $P[red] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.red,4)`. Similarly, $P[black] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.black,4)` and $P[green] = \frac{2}{38} =$ `r round(2/38,4)`, this is estimated as `r round(p.green,4)`

## Solution: Simulation Questions - Part 2

2. run the simulation with more trials, say 1 million?

```{r }
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.sims
```
Analytically $P[red] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.red,4)`. Similarly, $P[black] = \frac{18}{38} =$ `r round(18/38,4)`, this is estimated as `r round(p.black,4)` and $P[green] = \frac{2}{38} =$ `r round(2/38,4)`, this is estimated as `r round(p.green,4)`

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

