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-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.1V224H384V416H344V224H280V416H232V224H168V416H128V224zM256 64a32 32 0 1 1 0 64 32 32 0 1 1 0-64z"/></svg> EMSE 4571 / 6571: 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="M304 128a80 80 0 1 0 -160 0 80 80 0 1 0 160 0zM96 128a128 128 0 1 1 256 0A128 128 0 1 1 96 128zM49.3 464H398.7c-8.9-63.3-63.3-112-129-112H178.3c-65.7 0-120.1 48.7-129 112zM0 482.3C0 383.8 79.8 304 178.3 304h91.4C368.2 304 448 383.8 448 482.3c0 16.4-13.3 29.7-29.7 29.7H29.7C13.3 512 0 498.7 0 482.3z"/></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 24c0-13.3-10.7-24-24-24s-24 10.7-24 24V64H64C28.7 64 0 92.7 0 128v16 48V448c0 35.3 28.7 64 64 64H384c35.3 0 64-28.7 64-64V192 144 128c0-35.3-28.7-64-64-64H344V24c0-13.3-10.7-24-24-24s-24 10.7-24 24V64H152V24zM48 192H400V448c0 8.8-7.2 16-16 16H64c-8.8 0-16-7.2-16-16V192z"/></svg> April 25, 2024

]

---

# 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] "heads" "tails" "heads" "tails" "heads" "tails"
```

Probability of getting "heads":

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

```
#> [1] 0.5002
```

---

# 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] "heads" "tails" "heads" "heads" "tails" "tails"
```

Probability of getting "heads":

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

```
#> [1] 0.3915
```

---

# 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     2     6
#> 2     2     6     2
#> 3     6     6     5
#> 4     5     6     5
#> 5     4     3     5
#> 6     5     2     4
```

]

---

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

_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] 10  3  2
```

---

<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'
```

---

## Using `map` to take draws _without_ replacement

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

Make a function that returns 3 cards in a data frame:

```r
three_card_sum <- function(deck) {
  cards <- sample(x = deck, size = 3, replace = FALSE)
  return(sum(cards))
}

three_card_sum(deck)
```

```
#> [1] 26
```

]

Repeat the 3-card draw _N_ times:

```r
N <- 10000

sums <- map_int(1:N, function(x) three_card_sum(deck))

count <- sum(sums == 13)

count / N # Compute the probability
```

```
#> [1] 0.0359
```

]

---

## Using `for` loop to take draws _without_ replacement

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 <- 10000
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.0354
```

---

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

For face cards, assign numbers like this:

- 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[Intermission]

---

# Week 13: .fancy[Monte Carlo Methods]

### 1. Monte Carlo Simulation

### BREAK

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

---

# .center[Discrete vs. continuous random numbers]

### .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] "heads" "heads" "tails" "heads" "heads"
```
]

### .center[Continuous]

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

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

sample_continuous
```

```
#> [1] 0.4621284 0.7162149 0.8342098 0.4141555 0.4179559
```

]

---

# 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  6.32 23.3  TRUE      
#> 2  7.04  5.90 TRUE      
#> 3  5.63 24.1  TRUE      
#> 4  7.58 60.9  FALSE     
#> 5  7.43 46.3  TRUE      
#> 6  4.36 22.1  FALSE
```

]

]

---

**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.58527
```

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

```
#> [1] 149.8291
```
]

]

---

### How did we do?

Simulated area under curve:

```r
area_under_curve
```

```
#> [1] 149.8291
```

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.332
```

---

# 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 09, 3:00pm-5:00pm](https://my.gwu.edu/mod/exam_schedules/?x=E)