Monte Carlo Methods

]

# Week 13: .fancy[Monte Carlo Methods]

### <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M243.4 2.6l-224 96c-14 6-21.8 21-18.7 35.8S16.8 160 32 160v8c0 13.3 10.7 24 24 24H456c13.3 0 24-10.7 24-24v-8c15.2 0 28.3-10.7 31.3-25.6s-4.8-29.9-18.7-35.8l-224-96c-8.1-3.4-17.2-3.4-25.2 0zM128 224H64V420.3c-.6 .3-1.2 .7-1.8 1.1l-48 32c-11.7 7.8-17 22.4-12.9 35.9S17.9 512 32 512H480c14.1 0 26.5-9.2 30.6-22.7s-1.1-28.1-12.9-35.9l-48-32c-.6-.4-1.2-.7-1.8-1.1V224H384V416H344V224H280V416H232V224H168V416H128V224zm128-96c-17.7 0-32-14.3-32-32s14.3-32 32-32s32 14.3 32 32s-14.3 32-32 32z"/></svg> EMSE 4571: Intro to Programming for Analytics
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M272 304h-96C78.8 304 0 382.8 0 480c0 17.67 14.33 32 32 32h384c17.67 0 32-14.33 32-32C448 382.8 369.2 304 272 304zM48.99 464C56.89 400.9 110.8 352 176 352h96c65.16 0 119.1 48.95 127 112H48.99zM224 256c70.69 0 128-57.31 128-128c0-70.69-57.31-128-128-128S96 57.31 96 128C96 198.7 153.3 256 224 256zM224 48c44.11 0 80 35.89 80 80c0 44.11-35.89 80-80 80S144 172.1 144 128C144 83.89 179.9 48 224 48z"/></svg> John Paul Helveston
### <svg aria-hidden="true" role="img" viewBox="0 0 448 512" style="height:1em;width:0.88em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:white;overflow:visible;position:relative;"><path d="M152 64H296V24C296 10.75 306.7 0 320 0C333.3 0 344 10.75 344 24V64H384C419.3 64 448 92.65 448 128V448C448 483.3 419.3 512 384 512H64C28.65 512 0 483.3 0 448V128C0 92.65 28.65 64 64 64H104V24C104 10.75 114.7 0 128 0C141.3 0 152 10.75 152 24V64zM48 448C48 456.8 55.16 464 64 464H384C392.8 464 400 456.8 400 448V192H48V448z"/></svg> April 27, 2023

]

---

# Quiz 7

## Write your name on the quiz!

## Rules:

- Work alone; no outside help of any kind is allowed.
- No calculators, no notes, no books, no computers, no phones.

]

]

---

# Monte Carlo, Monaco

]]

]]

---

# "Monte Carlo" is associated with 3 things

### .center[Gambling]

]

### .center[Racing]

]

### .center[Simulation]

]

---

# Week 13: .fancy[Monte Carlo Methods]

### 1. Monte Carlo Simulation

### BREAK

### 2. Monte Carlo Integration

---

# Week 13: .fancy[Monte Carlo Methods]

### 1. .orange[Monte Carlo Simulation]

### BREAK

### 2. Monte Carlo Integration

---

# Monte Carlo Simulation: _Computing Probability_

### General process:

- Run a series of trials.
- In each trial, simulate an event (e.g. a coin toss, a dice roll, etc.).
- Count the number of "successful" trials
<br><br>

### `$\frac{\text{# Successful Trials}}{\text{# Total Trials}}$` = Observed Odds `$\simeq$` Expected Odds

### **Law of large numbers**:<br>As _N_ increases, Observed Odds <svg aria-hidden="true" role="img" viewBox="0 0 512 512" style="height:1em;width:1em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:currentColor;overflow:visible;position:relative;"><path d="M470.6 278.6c12.5-12.5 12.5-32.8 0-45.3l-160-160c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3L402.7 256 265.4 393.4c-12.5 12.5-12.5 32.8 0 45.3s32.8 12.5 45.3 0l160-160zm-352 160l160-160c12.5-12.5 12.5-32.8 0-45.3l-160-160c-12.5-12.5-32.8-12.5-45.3 0s-12.5 32.8 0 45.3L210.7 256 73.4 393.4c-12.5 12.5-12.5 32.8 0 45.3s32.8 12.5 45.3 0z"/></svg> Expected Odds

---

# How would you measure if a coin is "fair"?

Run a series of trials and record outcome: "heads" or "tails"

```r
coin <- c("heads", "tails")
N <- 10000
*tosses <- sample(x = coin, size = N, replace = TRUE)
head(tosses) # Preview first few tosses
```

```
#> [1] "tails" "tails" "tails" "heads" "tails" "heads"
```

Probability of getting "heads":

```r
sum(tosses == "heads") / N
```

```
#> [1] 0.5018
```

---

# Tossing an unfair coin

Set the `prob` argument to a 40-60 coin

```r
coin <- c("heads", "tails")
N <- 10000
*tosses <- sample(x = coin, size = N, replace = TRUE, prob = c(0.4, 0.6))
head(tosses) # Preview first few tosses
```

```
#> [1] "tails" "heads" "tails" "heads" "heads" "tails"
```

Probability of getting "heads":

```r
sum(tosses == "heads") / N
```

```
#> [1] 0.4026
```

---

# A more complex simulation: _dice rolling_

### What is the probability of rolling a 6-sided dice 3 times<br>and getting the sequence 1, 3, 5?

```r
library(tidyverse)
dice <- c(1, 2, 3, 4, 5, 6)
N <- 10000
rolls <- tibble(
  roll1 = sample(x = dice, size = N, replace = T),
  roll2 = sample(x = dice, size = N, replace = T),
  roll3 = sample(x = dice, size = N, replace = T)
)
```

]]

```r
head(rolls)
```

```
#> # A tibble: 6 × 3
#>   roll1 roll2 roll3
#>   <dbl> <dbl> <dbl>
#> 1     4     5     2
#> 2     4     1     1
#> 3     1     5     4
#> 4     2     6     6
#> 5     1     3     5
#> 6     6     2     6
```

]

---

# A more complex simulation: _dice rolling_

Simulated probability of getting sequence 1, 3, 5:

```r
successes <- rolls %>%
  filter(roll1 == 1 & roll2 == 3 & roll3 == 5)

nrow(successes) / N
```

```
#> [1] 0.0049
```

_Actual_ probability of getting sequence 1, 3, 5:

```r
(1/6)^3
```

```
#> [1] 0.00462963
```

---

## Your Turn: Coins & Dice

Using the `sample()` function, conduct a monte carlo simulation to estimate the answers to these questions:

- If you flipped a coin 3 times in a row, what is the probability that you'll get three "tails" in a row?
- If you rolled 2 dice, what is the probability that you'll get "snake-eyes" (two 1's)?
- If you rolled 2 dice, what is the probability that you'll get an outcome that sums to 8?

]

---

Sometimes events cannot be independently simulated

<br>

What are the odds that 3 cards drawn from a 52-card deck will sum to 13?

- Aces = 1
- Jack = 10
- Queen = 10
- King = 10

---

Sometimes events cannot be independently simulated

<br>

```r
cards <- c(seq(1, 10), 10, 10, 10)
deck <- rep(cards, 4) # Rep because there are 4 suits
length(deck)
```

```
#> [1] 52
```

Draw 3 cards from the deck _without replacement_:

```r
*draw <- sample(x = deck, size = 3, replace = FALSE)
draw
```

```
#> [1]  9  4 10
```

---

<br>

**Note**: You can't draw more than 52 cards _without replacement_:

```r
draw <- sample(x = deck, size = 53, replace = FALSE)
```

```
#> Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'
```

---

What are the odds that 3 cards drawn from a 52-card deck will sum to 13?

<br>

Repeat the 3-card draw _N_ times:

```r
N <- 100000
count <- 0

for (i in 1:N) {
*  draw <- sample(x = deck, size = 3, replace = FALSE)
   if (sum(draw) == 13) {
      count <- count + 1
   }
}

count / N # Compute the probability
```

```
#> [1] 0.03662
```

---

## Your Turn: Cards

Use the `sample()` function and a monte carlo simulation to estimate the answers to these questions:

- What are the odds that four cards drawn from a 52-card deck will have the same suit?
- What are the odds that five cards drawn from a 52-card deck will sum to a prime number?

- Aces = 1
- Jack = 10
- Queen = 10
- King = 10

]

**Hint**: use `isPrime()` to help:

```r
isPrime <- function(n) {
    if (n == 2) { return(TRUE) }
    for (i in seq(2, n-1)) {
        if (n %% i == 0) {
            return(FALSE)
        }
    }
    return(TRUE)
}
```

]

---

# .fancy[Break]

---

# Week 13: .fancy[Monte Carlo Methods]

### 1. Monte Carlo Simulation

### BREAK

### 2. .orange[Monte Carlo Integration]

---

### .center[Discrete]

`sample()`<br>Takes random samples from vector `x`

```r
sample_discrete <- sample(
  x       = c("heads", "tails"),
  size    = 5,
  replace = TRUE
)

sample_discrete
```

```
#> [1] "tails" "heads" "heads" "heads" "tails"
```
]

### .center[Continuous]

`runif()`<br>Takes random samples between bounds

```r
sample_continuous <- runif(
  n   = 5,
  min = 0,
  max = 1
)

sample_continuous
```

```
#> [1] 0.7536088 0.2830033 0.6428447 0.4526490 0.7055930
```

]

---

# Monte Carlo Integration

Integration = compute the area "under the curve"

Find the area of `$y = x^2$` between `$4 < x < 8$`

]]

`$\frac{\text{Area Under Curve}}{\text{Area of Rectangle}} = \frac{\text{# Points Under Curve}}{\text{# Total Points}}$`
<br>

]

---

`$\frac{\text{Area Under Curve}}{\text{Area of Rectangle}} = \frac{\text{# Points Under Curve}}{\text{# Total Points}}$`

`$\text{Area Under Curve} = \text{Area of Rectangle} \left( \frac{\text{# Points Under Curve}}{\text{# Total Points}} \right)$`

---

**Step 1**: Compute area of rectangle

```r
area_rectangle <- (8 - 4) * (8^2 - 0)
area_rectangle
```

```
#> [1] 256
```
]

]

---

**Step 2**: Simulate points

```r
N <- 100000
points <- tibble(
*   x = runif(N, min = 4, max = 8),
*   y = runif(N, min = 0, max = 8^2)) %>%
    mutate(belowCurve = y < x^2)

head(points)
```

```
#> # A tibble: 6 × 3
#>       x     y belowCurve
#>   <dbl> <dbl> <lgl>     
#> 1  7.59 20.6  TRUE      
#> 2  4.51 25.7  FALSE     
#> 3  7.46  4.37 TRUE      
#> 4  4.73 56.4  FALSE     
#> 5  6.43 27.1  TRUE      
#> 6  7.35 48.0  TRUE
```

]

]

---

**Step 3**: Compute area under curve

```r
N <- 100000
points <- tibble(
    x = runif(N, min = 4, max = 8),
    y = runif(N, min = 0, max = 8^2)) %>%
    mutate(belowCurve = y < x^2)

*points_ratio <- sum(points$belowCurve) / N
points_ratio
```

```
#> [1] 0.58365
```

```r
*area_under_curve <- area_rectangle * points_ratio
area_under_curve
```

```
#> [1] 149.4144
```
]

]

---

### How did we do?

Simulated area under curve:

```r
area_under_curve
```

```
#> [1] 149.4144
```

Actual area under curve:

`$\int_{4}^{8} x^2 \mathrm{dx} = \left ( \frac{x^3}{3} \right ) \Big|_4^8 = \frac{8^3}{3} - \frac{4^3}{3} = 149.33\bar{3}$`

% Error:

```r
true_area <- ((8^3 / 3) - (4^3 / 3))
100*((area_under_curve - true_area) / true_area)
```

```
#> [1] 0.05428571
```

---

# Monte Carlo `$\pi$`

]

Area of a circle:

`$\quad\quad A_{circle} = \pi r^2$`

Area of square containing circle:

`$\quad\quad A_{square} = 4r^2$`

]

---

# Monte Carlo `$\pi$`

]

Area of a circle:

`$\quad\quad A_{circle} = \pi r^2$`

Area of square containing circle:

`$\quad\quad A_{square} = 4r^2$`

Ratio of areas = `$\pi / 4$`:

`$\quad\quad \frac{A_{circle}}{A_{square}} = \dfrac{\pi r^2}{4r^2} = \dfrac{\pi}{4}$`

]

---

# Monte Carlo `$\pi$`

]

Area of a circle:

`$\quad\quad A_{circle} = \pi r^2$`

Area of square containing circle:

`$\quad\quad A_{square} = 4r^2$`

Ratio of areas = `$\pi / 4$`:

`$\quad\quad \frac{A_{circle}}{A_{square}} = \dfrac{\pi r^2}{4r^2} = \dfrac{\pi}{4}$`

`$\pi = 4 \left( \frac{A_{circle}}{A_{square}} \right)$`

]

---

## Your Turn: Estimate `$\pi$`

`$\pi = 4 \left( \frac{A_{circle}}{A_{square}} \right)$`

]
.rightcol70[.font90[

1. Create a tibble with variables `x` and `y` that each contain 10,000 random points between -1 and 1, representing the (x, y) coordinates to a random point inside a square of side length 2 centered at `(x, y) = (0, 0)`. **Hint**: use `runif()`

2. Create a new column, `radius`, that is equal to the distance to each `(x, y)` point from the center of the square.

3. Create a new column, `pointInCircle`, that is `TRUE` if the point lies _within_ the circle inscribed in the square, and `FALSE` otherwise.

4. Create the scatterplot on the left (don't worry about the precise colors, dimensions, etc.).

5. Estimate `$\pi$` by multiplying 4 times the ratio of points inside the circle to the total number of points

]]

---

---

## Your Turn: Monte Hall Problem

1. You choose door 1, 2, or 3
2. One door is removed
3. Should you swap doors?

]

In this simulation, the prize is always behind door #1:

- If you choose door #1, you must KEEP it to win.
- If you choose door #2 or #3, you must SWAP to win.

1) Create the tibble, `choices`, with two variables:
- `door` contains the first door chosen (`1`, `2`, or `3`)
- `swap` contains a logical (`TRUE` or `FALSE`) for whether the contestant swaps doors. **Hint**: use `sample()`

2) Create a new tibble, `wins`, which contains only the rows from `choices` that resulted in a win.

3) Compute the percentage of times the contestant won after swapping doors.

]]

---

# Reminders

### 1) Please fill the GW course feedback (see slack announcement)

### 2) Final is [Thursday, May 11, 3:00pm-5:00pm](https://my.gwu.edu/mod/exam_schedules/?x=E)