Statistical Simulation

Exploring What-If Scenarios

Chapter 4

Why Simulate?

The Central Question

How might our plans, methods, model \(\ldots\) fail?

Simulation lets us explore scenarios we can’t observe directly.

Generate data from known processes → see how methods behave.

Uses of Simulation

  1. Compare methods: Which estimator performs better?
  1. Assess robustness: What happens with outliers? Small samples?
  1. Understand systems: How does randomness propagate?
  1. Quantify uncertainty: How precise are our estimates?

The EDA Connection

Simulation extends EDA beyond the data we have.

“What would happen if \(\dots\)?”

  • \(\dots\) the sample size were larger?
  • \(\dots\) the distribution had heavier tails?
  • \(\dots\) there were outliers?

Simulation lets us explore failure modes before they occur in practice.

Random Number Generation

R’s Distribution Functions

For each distribution, R provides four functions:

Prefix Function Example
d density/mass dnorm(0) → 0.399
p cumulative probability pnorm(1.96) → 0.975
q quantile (inverse CDF) qnorm(0.975) → 1.96
r random generation rnorm(10) → 10 values

Continuous Distributions (selected)

Table 1
Function Distribution
[d,p,q,r]norm Normal
[d,p,q,r]unif Uniform
[d,p,q,r]exp Exponential
[d,p,q,r]gamma Gamma
[d,p,q,r]beta Beta
[d,p,q,r]t Student t
[d,p,q,r]chisq Chi-squared

See ?Distributions for a complete list.

Discrete Distributions (selected)

Table 2
Function Distribution
[d,p,q,r]binom Binomial
[d,p,q,r]pois Poisson
[d,p,q,r]geom Geometric
[d,p,q,r]nbinom Negative Binomial

The Poisson is the “law of rare events”—the limit of Binomial\((n, p)\) as \(n \to \infty\) and \(p \to 0\) with \(np \to \lambda\).

Generating Random Samples

Code 
# 100 standard normal values
x <- rnorm(n = 100, mean = 0, sd = 1)

# 100 uniform values on [0, 1]
u <- runif(n = 100, min = 0, max = 1)

# 100 Poisson values with mean 5
counts <- rpois(n = 100, lambda = 5)

Key insight: We specify the distribution; R generates samples.

Comparing Estimators

Mean vs. Median

Both estimate the center of a symmetric distribution.

Mean: \(\hat{\mu} = \frac{1}{n}\sum_{i=1}^n X_i\)

Median: \(\hat{m}\) = middle value (or average of two middle values)

Which is better? It depends on the distribution.

The Efficiency Question

For normal data, theory tells us:

\[\frac{\text{Var}(\hat{m})}{\text{Var}(\hat{\mu})} \approx \frac{\pi}{2} \approx 1.57\]

The median has 57% higher variance than the mean.

Can we verify this by simulation?

Simulation Design

Code 
# Parameters
n <- 30      # sample size
R <- 10000   # number of replications

# Storage
means <- numeric(R)
medians <- numeric(R)

# Simulation loop
for (r in 1:R) {
  x <- rnorm(n)
  means[r] <- mean(x)
  medians[r] <- median(x)
}

# Compare variances
var(medians) / var(means)

Simulation Results

Figure 1: Sampling distributions of mean and median (n=30, R=10,000)

Verifying Theory

Table 3
Estimator Variance Std Error
Mean 0.0338 0.1837
Median 0.0512 0.2263

Variance ratio: 1.517 (theory: \(\pi/2 \approx\) 1.571)

Simulation confirms the theoretical result.

Why This Matters

We verified a known result. But simulation shines when:

  • Theory is intractable
  • Assumptions are violated
  • We need to compare many methods

EDA application: “How would my estimator behave if the data were slightly different?”

Monte Carlo Methods

Estimating π

Classic example: estimate \(\pi\) by random sampling.

Idea: A quarter circle of radius 1 has area \(\pi/4\).

  1. Generate random points \((X, Y)\) uniform on \([0,1]^2\)
  2. Count points inside the quarter circle: \(X^2 + Y^2 < 1\)
  3. Proportion inside \(\approx \pi/4\)

Visual Intuition

Figure 2: Monte Carlo estimation of π

The Monte Carlo Principle

To estimate \(E[g(X)]\):

\[\hat{\theta} = \frac{1}{n}\sum_{i=1}^n g(X_i)\]

where \(X_1, \ldots, X_n\) are random samples from the distribution of \(X\).

Law of Large Numbers: \(\hat{\theta} \to E[g(X)]\) as \(n \to \infty\)

Convergence Rate

Figure 3: Monte Carlo estimate converges to true value

Standard error decreases as \(1/\sqrt{n}\).

Rare Events

The Problem

What if we want to estimate \(P(Z > 6)\) where \(Z \sim N(0,1)\)?

Code 
# True probability
pnorm(6, lower.tail = FALSE)
[1] 9.865876e-10

About 1 in a billion. Naive simulation won’t work.

Naive Approach Fails

Code 
# Generate 1 million standard normals
z <- rnorm(1e6)

# Count how many exceed 6
sum(z > 6)
[1] 0

We’d need billions of samples to see even one event.

Importance Sampling

Key idea: Sample from a different distribution, then reweight.

Instead of sampling from \(p(x)\), sample from \(q(x)\) where rare events are common.

\[E_p[g(X)] = E_q\left[g(X) \cdot \frac{p(X)}{q(X)}\right]\]

The ratio \(w(x) = p(x)/q(x)\) is the importance weight.

Importance Sampling for \(P(Z > 6)\)

Strategy: Shift the normal distribution to center it near 6.

Let \(q(x) = \phi(x - 6)\) (normal centered at 6).

Code 
# Sample from shifted normal
n <- 10000
x <- rnorm(n, mean = 6, sd = 1)

# Importance weights: p(x) / q(x)
log_weights <- dnorm(x, log = TRUE) - dnorm(x, mean = 6, log = TRUE)
weights <- exp(log_weights)

# Estimate P(Z > 6)
estimate <- mean((x > 6) * weights)
estimate
[1] 9.91751e-10

Comparing Estimates

Table 4
Method Estimate
True value 9.87e-10
Importance sampling 9.92e-10
Naive (n = 10^6) 0

Importance sampling makes the impossible feasible.

Bootstrap

Quantifying Uncertainty

How precise is our estimate?

Classical approach: derive standard error formula.

Bootstrap approach: Resample from the data itself.

The Bootstrap Idea

  1. Treat the sample as a proxy for the population
  2. Draw with replacement from the sample
  3. Compute statistic on each resample
  4. Use variability across resamples to estimate uncertainty

Bootstrap in Action

Code 
# Original sample
x <- c(12, 15, 18, 22, 25, 28, 31, 35, 42, 55)

# Bootstrap
B <- 10000
boot_means <- numeric(B)

for (b in 1:B) {
  x_star <- sample(x, replace = TRUE)
  boot_means[b] <- mean(x_star)
}

# 95% confidence interval
quantile(boot_means, c(0.025, 0.975))

Bootstrap Example

Figure 4: Bootstrap distribution of the sample mean

Normal theory 95% CI: [18.9, 37.7] — similar result.

Bootstrap for Any Statistic

The bootstrap works for almost any statistic:

  • Correlation coefficients
  • Regression coefficients
  • Median, trimmed mean
  • Ratios, differences

EDA application: “How much would this pattern change with different data?”

Connections to ML

Simulation in the ML Pipeline

Simulation connects to machine learning:

Cross-validation: Repeatedly resample to estimate generalization error.

Stochastic gradient descent: Random sampling of training batches.

Bayesian ML: MCMC and variational inference sample from posteriors.

Simulation as Model Checking

After fitting a model:

  1. Simulate data from the fitted model
  2. Compare simulated data to real data
  3. Discrepancies reveal model failures

This is the posterior predictive check—simulation for model criticism.

Summary

Chapter 4: Key Takeaways

  1. Simulation explores what-if scenarios before they occur

  2. R provides d/p/q/r functions for common distributions

  3. Monte Carlo methods estimate quantities via random sampling

  4. Importance sampling handles rare events

  5. Bootstrap quantifies uncertainty without formulas

Key Formulas

Concept Formula
Monte Carlo estimate \(\hat{\theta} = \frac{1}{n}\sum_{i=1}^n g(X_i)\)
Importance weight \(w(x) = p(x) / q(x)\)
Variance ratio (mean vs. median) \(\pi/2 \approx 1.57\)

The Simulation Workflow

  1. Specify the data-generating process
  2. Generate random samples
  3. Compute statistics of interest
  4. Repeat many times
  5. Summarize the distribution of results

Exercises

Team Exercise 1: Mean vs. Median Efficiency

Verify the variance ratio for normal data:

  1. Set \(n = 30\) and \(R = 5000\) replications.
  2. For each replication, generate a sample from rnorm(), compute mean and median.
  3. Calculate var(medians) / var(means). Compare to \(\pi/2 \approx 1.57\).
  4. What does this ratio tell us about the relative efficiency of the median?

Team Exercise 2: Heavy Tails and the Cauchy

Repeat Exercise 1, but use rcauchy() instead of rnorm().

  1. Calculate the sample variance of your \(R\) means. What do you observe?
  2. Construct a histogram of sample means. Does it look normal?
  3. What does this imply about the Central Limit Theorem?

Team Exercise 3: Bootstrap Confidence Interval

Use the bootstrap to estimate a 95% CI for the median of mtcars$mpg:

  1. Resample with replacement \(B = 2000\) times.
  2. Compute the median for each bootstrap sample.
  3. Use the 2.5th and 97.5th percentiles as your CI.
  4. How does the width compare to a bootstrap CI for the mean?

Team Exercise 4: Importance Sampling

Estimate \(P(Z > 4)\) using importance sampling. Compare to pnorm(4, lower.tail = FALSE).

Discussion Questions

  1. When is simulation more trustworthy than mathematical derivation?

  2. A colleague sets set.seed(42) and runs one simulation. Is this reproducible research?

  3. How many bootstrap replications are “enough”?

Resources