STAT 302 Statistical Computing

class: center, top, title-slide

.title[
# STAT 302 Statistical Computing
]
.subtitle[
## Lecture 6: Simulations
]
.author[
### Yikun Zhang (Winter 2024)
]

---

# Outline

1. Simulation Basics

2. Pseudo-Randomness and Seeds

3. Simulation Tools in R

4. Integral Approximation via Monte Carlo Methods

5. Distribution Estimation and Sampling via Monte Carlo Methods

6. Bootstrapping

7. Ad hoc Network (Final Project)

* Acknowledgement: Parts of the slides are modified from the course materials by Prof. Ryan Tibshirani, Prof. Yen-Chi Chen, Prof. Alexander Giessing, Prof. Deborah Nolan, and Prof. Noureddine El Karoui.

---
class: inverse

# Part 1: Simulation Basics

---

# What is Simulation?

- A simulation imitates the operation of real-world processes or systems with the use of a computer software.

- In most cases, simulation involves creating data by (pseudo-)random number generators and using that data to study a problem of interest.
  
  - Such simulation process is also known as the [Monte Carlo Methods](https://en.wikipedia.org/wiki/Monte_Carlo_method), whose name originates from the Monte Carlo Casino in Monaco.
  
--
  
Can you think of any examples from this or other classes when you have used simulations to answer a question?

---

# An Example of Simulation in Cosmology

<img src="./figures/illustris_z0.png" width="400"/>


Exterior view of the dark matter density distribution in the full simulation box at redshift zero for <a href="https://www.illustris-project.org/">Illustris Simulation Project</a>.

---

# Why Do We Need Simulations?

To answer some scientific questions, we often need data evidence and experimental results to support our findings.

- However, observational data are generally expensive to obtain.

- On the contrary, it is relatively at a lower cost to generate synthetic data and implement a computer programming to conduct experiments.

- In addition, simulations sometimes can be easier and more realistic than hand calculations.

---

# Concrete Motivations for Simulations

A well-designed simulation can provide a solid estimator to statistical quantities and their uncertainty measures that are analytically difficult to obtain the exact answer.

- What is the probability of flipping 10 fair coins and getting 9 heads and 1 tails?

- We can compute by hand as `${10 \choose 9} \left(\frac{1}{2}\right)^9 \left(1-\frac{1}{2}\right) \approx 0.0097656$`.
  
--
  
  - Nevertheless, the simulation-based approximation could be more convenient.

```r
set.seed(123)
coins = rbinom(n = 5000000, size = 10, prob = 0.5)
mean(coins == 9)
```

```
## [1] 0.0097274
```

---

# Concrete Motivations for Simulations

It also helps evaluate the performance of our proposed statistical method, such as its asymptotic normality and time efficiency.

Simulation can also lead to a more efficient hypothesis testing or other statistical procedures than the original one.

- Example: The Monte Carlo alternative to the [permutation test](https://en.wikipedia.org/wiki/Permutation_test).

---

# Concrete Motivations for Simulations
  
Simulation also provides a feasible way to generate random samples from a "partially known" probability distribution.

<img src="./figures/rejectionsampling.png" width="550"/>


 <a href="https://en.wikipedia.org/wiki/Rejection_sampling">Acceptance-Rejection Sampling</a>

---

# Concrete Motivations for Simulations

Does simulation only appear in Statistics? Can other fields benefit from the usage of simulation?

- Monte Carlo methods have been developed into a technique called [Monte Carlo tree search](https://en.wikipedia.org/wiki/Monte_Carlo_tree_search) that is useful for searching for the best move in a game.

- The Monte Carlo tree search has been successfully combined with deep neural networks to tackle the multiple board games, such as Go, Chess, Texas hold 'em, etc. See the introduction to [AlphaGo](https://www.deepmind.com/research/highlighted-research/alphago).

---

# Overview of Monte Carlo Simulations

A typical Monte Carlo simulation often follows a procedure as below.

1. Define a domain of possible inputs.

2. Generate inputs randomly from a probability distribution over the domain.

3. Perform a deterministic computation on the inputs.

4. Aggregate the results.

Do you know what the first known Monte Carlo simulation is?

<img src="./figures/Buffon_needle.png" width="300"/>


 <a href="https://en.wikipedia.org/wiki/Buffon%27s_needle_problem">Buffon's Needle Problem</a> for Estimating `$\pi$`

---

# An Example of Monte Carlo Simulation

```r
set.seed(123)
n = 400000
x = runif(n, min = 0, max = 1)
y = runif(n, min = 0, max = 1)
mean(x^2 + y^2 < 1) * 4
```

```
## [1] 3.14125
```

---
class: inverse

# Part 2: Pseudo-Randomness and Seeds

---

# Random Number Generation in R

Most of the simulations require us to generate some random observations from a probability distribution.

Recall from Lab 2 that we have the utility functions:

- `rnorm()`: generates the normally distributed random variables/observations.

- `rbinom()`: generates the random variables/observations from a binomial distribution.

We can replace "norm" or "binom" with the name of another distribution, and the same function applies, e.g., "t", "unif", "exp", "gamma", "chisq", "pois", etc.

---

# How Does R Generate Random Numbers?

There are no computer software that can generate completely random numbers, neither did R.

Instead, R uses a **pseudo-random number generator**:

- It starts with a **seed** and an **algorithm** (i.e., a function).

- The **seed** is input into the algorithm, and a (pseudo-random) number is returned. This number is then plugged into the algorithm, and the next (pseudo-random) number is created.

- The **algorithm** guarantees that the output numbers behave like random values.

---

# Linear Congruential Generator

One example of such an algorithm is the **Linear Congruential Generator**, which uses modular arithmetic to generate "random" numbers:

- Based on two parameters `$m$` and `$b$` and the initial value `$x_0$` (i.e., the seed), the first "random number" `$x_1$` is generated as follows:
`$$x_1 = b \cdot x_0 \mod m.$$`
- The subsequent numbers are generated recursively,
`$$x_{n+1} = b \cdot x_n \mod m$$`
for `$n=1,2,...$`.

---

# Linear Congruential Generator (Small `$m$`)

```r
LinCongGen = function(n, seed, m = 64, b = 3){
  pseudo_num = numeric(n+1)
  pseudo_num[1] = seed
  for(i in 1:n){
    pseudo_num[i+1] = (pseudo_num[i]*b) %% m
  }
  return(pseudo_num[2:(n+1)])
}
n = 1000
seed = 17
rand_num = LinCongGen(n, seed, m = 64, b = 3)
plot(c(seed, rand_num[1:(n-1)]), rand_num, xlab = "Input Numbers", ylab = "Output Numbers")
```

---

# Linear Congruential Generator (Large `$m$`)

```r
n = 1000
seed = 17
rand_num = LinCongGen(n, seed, m = 2^32, b = 3)
plot(c(seed, rand_num[1:(n-1)]), rand_num, xlab = "Input Numbers", ylab = "Output Numbers")
```

The output numbers are still not "random" enough.

---

# Linear Congruential Generator (Large `$m$`, Large `$b$`)

```r
n = 1000
seed = 17
rand_num = LinCongGen(n, seed, m = 2^32, b = 69069)
plot(c(seed, rand_num[1:(n-1)]), rand_num, xlab = "Input Numbers", ylab = "Output Numbers")
```

---

# Pseudo-Random Number Generator in R

- Studying algorithms of generating more "random" numbers is an interesting research area in its own right.

- The default algorithm in R (and in nearly all software languages) is called the "Mersenne-Twister" algorithm.

- Type `?Random` in our R console to read more about this (and to read how to change the algorithm used for the pseudo-random number generation, which we should never really have to do that).

---

# An Important Role of the "Seed"

One of the biggest advantages of pseudo-random number generators is that we can have some controls over the simulation results by setting the **seed** beforehand.

- The seed is just an integer, and can be set with `set.seed()`.

- It is very important to set our seed in the **proper** location of our code to make our simulation results _valid_ and _reproducible_.

---

# Demonstrations of `set.seed()`

If we don't set the seed, it is not surprising to see different draws each time we call any random number generating functions.

```r
mean(rnorm(6))
```

```
## [1] 0.2484921
```

```r
mean(rnorm(6))
```

```
## [1] -0.6256406
```

```r
set.seed(123)
mean(rnorm(6))
```

```
## [1] 0.4471527
```

```r
set.seed((123))
mean(rnorm(6))
```

```
## [1] 0.4471527
```

---

# Demonstrations of `set.seed()`

Each time the seed is set, the same sequence follows (indefinitely).

```r
set.seed(123)
rnorm(3); rnorm(2); rnorm(1)
```

```
## [1] -0.5604756 -0.2301775  1.5587083
```

```
## [1] 0.07050839 0.12928774
```

```
## [1] 1.715065
```

```r
set.seed(123)
rnorm(3); rnorm(2); rnorm(1)
```

```
## [1] -0.5604756 -0.2301775  1.5587083
```

```
## [1] 0.07050839 0.12928774
```

```
## [1] 1.715065
```

---
class: inverse

# Part 3: Simulation Tools in R

---

# `sample()` Function

The `sample()` function takes a sample of the pre-specified size from the elements of `x` using either with or without replacement.

```r
sample(letters, size = 10, replace = FALSE, prob = NULL)
```

```
##  [1] "v" "y" "e" "s" "i" "c" "h" "g" "j" "z"
```

We can make the sampling results reproducible by setting the seed.

```r
set.seed(123)
sample(letters, size = 10, replace = FALSE, prob = NULL)
```

```
##  [1] "o" "s" "n" "c" "j" "r" "k" "e" "x" "y"
```

---

# Betting Red in Roulette

What is the chance that a spin of the Roulette wheel lands on red?

- There are 38 pockets in the Roulette wheel:
  
  - 18 blacks, 18 reds, and 2 greens.

- We can estimate the chance by simulating the Roulette wheel many, many, many times via the `sample()` function in R:
`$$\frac{\# \text{ reds}}{\# \text{ draws}} \approx \text{Chance of reds}.$$`

---

# Betting Red in Roulette

What is the chance that a spin of the Roulette wheel lands on red?

```r
rou_vec = c(rep(c("black", "red"), each = 18), rep("green", 2))

set.seed(123)
samp_res = sample(rou_vec, size = 10000, replace = TRUE)
mean(samp_res == "red")
```

```
## [1] 0.4759
```

```r
cat("The estimation error is ", abs(mean(samp_res == "red") - 18/38), ".", sep = "")
```

```
## The estimation error is 0.002215789.
```

---

# Betting Red in Roulette

As the number of draws increases, the proportions of reds becomes closer to the true probability `$\frac{18}{38}$`.

```r
set.seed(123)
samp_res = sample(rou_vec, size = 1000000, replace = TRUE)
cat("The estimation error is ", abs(mean(samp_res == "red") - 18/38), ".", sep = "")
```

```
## The estimation error is 0.0005752105.
```

This reveals the principle of **Law of Large Numbers** in Statistics:
`$$\frac{1}{n} \sum_{i=1}^n X_i \to p \quad \text{ in probability } \quad \text{or even} \quad \text{ almost surely},$$`
where `$p$` is the probability of lands on reds in the Roulette wheel, and 
$$
X_i = 
`\begin{cases}
1 & \text{ when landing on reds for the i-th draw},\\
0 & \text{ when landing on other colors for the i-th draw}.\\
\end{cases}`
$$

---

# Betting Red in Roulette

By increasing the sample size, the precision of the estimated probability of the Roulette wheel landing on reds can be improved.

- However, for any fixed sample size `$n$`, what are the typical deviations (or standard errors) from the average of these sample proportions of reds?

- To quantify these deviations (or standard errors), we need to repeat the previous sampling procedure `$B$` times, where `$B$` is sufficiently large.
  
  - In R, such repetitions can be handled by the `replicate()` function.

---

# `replicate()` Function

Let's fix the sample size `$n$` to be 10000 and repeat the entire procedure of estimating the probability of the Roulette wheel landing on reds `$B=3000$` times.

```r
set.seed(123)   ## We should set the seed here!!
rep_means = replicate(3000, 
                      mean(sample(rou_vec, size = 10000, replace = TRUE) == "red"))
```

Note: The `replicate` function is a wrapper for the common use of `sapply()` for repeated evaluation of an expression.

```r
set.seed(123)   ## We should set the seed here!!
rep_means2 = replicate(3000, {
  samp_res = sample(rou_vec, size = 10000, replace = TRUE)
  mean(samp_res == "red")
})
all.equal(rep_means, rep_means2)
```

```
## [1] TRUE
```

---

# Distribution of the Sample Proportions

```r
hist(rep_means, breaks = 30, freq = FALSE, xlab = "Sample Proportion of Reds", 
     main = "Histogram of the Sample Proportions")
```

---

# Rate of Convergence of the Sample Proportions

How quickly do the sample proportions approach the true probability `$\frac{18}{38}$` as the sample size `$n$` increases?

Note: The red curves plot the trends `$\pm \sqrt{\frac{p(1-p)}{n}}$`.

---

# Distribution of the Sample Proportions

- The distribution of the sample proportions of reds looks roughly like the normal distribution with

- Center: `mean`(sample proportions) `$\approx \frac{18}{38}$`.
  
  - Spread: `sd`(sample proportions) `$\approx \sqrt{\frac{p(1-p)}{n}}$`.
  
---

# Distribution of the Sample Proportions

This reveals the **Central Limit Theorem** in Statistics:
`$$\frac{\sqrt{n} (\bar{X}_n - p) }{\sqrt{p(1-p)}} \to N(0,1) \quad \text{ in distribution},$$`
where `$\bar{X}_n = \frac{1}{n} \sum_{i=1}^n X_i$`.

---