5.2 Simulation & Monte Carlo Methods

Why this week matters¶

Week 9 focused on vectorized numerical computing with NumPy. Week 10 uses that machinery to study one of the most important tools in computational statistics:

simulation

The approved STA 556 schedule identifies Week 10 as covering:

random-number generation;
inverse-transform sampling;
Monte Carlo methods.

Simulation is fundamental because many statistical quantities are difficult or impossible to derive analytically, but can be approximated computationally.

Typical uses include:

approximating probabilities;
estimating expectations;
evaluating estimators;
exploring sampling distributions;
approximating integrals;
propagating uncertainty;
validating statistical procedures.

1. What is simulation?¶

Simulation means generating artificial observations from a specified probabilistic model.

For example:

import numpy as np

rng = np.random.default_rng(556)

x = rng.normal(
    loc=0,
    scale=1,
    size=1000
)

Here we simulate:

X₁, ..., X₁₀₀₀ ~ Normal(0,1)

The computer generates values that behave like draws from the model.

2. Pseudorandom numbers¶

Computers do not usually generate truly random numbers algorithmically.

Instead, they generate pseudorandom sequences.

A pseudorandom generator:

seed
  ↓
deterministic algorithm
  ↓
sequence that behaves statistically like random draws

This is valuable because the same seed can reproduce the same simulation.

3. Reproducible random-number generation¶

Prefer NumPy’s generator API:

rng = np.random.default_rng(556)

Then:

rng.normal(size=5)

If the same seed is used again:

rng = np.random.default_rng(556)

the same sequence is generated.

This makes simulation reproducible.

4. Why `default_rng()`?¶

Modern NumPy recommends generator objects rather than relying on global random state.

A generator makes randomness an explicit dependency.

Instead of:

np.random.seed(556)
np.random.normal(size=10)

prefer:

rng = np.random.default_rng(556)
rng.normal(size=10)

The generator can also be passed into functions.

5. Randomness as a function dependency¶

Consider:

def simulate_scores(
    n,
    mean,
    sd,
    rng
):
    return rng.normal(
        mean,
        sd,
        size=n
    )

Then:

rng = np.random.default_rng(556)

scores = simulate_scores(
    n=100,
    mean=80,
    sd=10,
    rng=rng
)

The source of randomness is explicit.

This makes the function easier to test and reproduce.

6. Common distributions¶

NumPy generators support many useful distributions.

Normal:

rng.normal(
    loc=0,
    scale=1,
    size=1000
)

Uniform:

rng.uniform(
    low=0,
    high=1,
    size=1000
)

Binomial:

rng.binomial(
    n=10,
    p=0.4,
    size=1000
)

Poisson:

rng.poisson(
    lam=3,
    size=1000
)

Choice from categories:

rng.choice(
    ["A", "B", "C"],
    size=100
)

7. Verify a simulation¶

If:

x = rng.normal(
    10,
    2,
    size=100000
)

then check:

x.mean()
x.std()

We expect approximately:

mean ≈ 10
sd ≈ 2

Simulation output should be checked against theoretical expectations.

8. Simulation variability¶

If we simulate only:

n = 10

the sample mean may differ considerably from the population mean.

With:

n = 100000

the sample mean will usually be closer.

Simulation results themselves are random.

This leads to a central Monte Carlo idea:

Approximation error decreases as the number of simulations increases.

9. Monte Carlo methods¶

A Monte Carlo method uses repeated random simulation to approximate a numerical or probabilistic quantity.

Generic workflow:

Define target quantity
       ↓
Generate random samples
       ↓
Compute statistic
       ↓
Repeat many times
       ↓
Average / summarize
       ↓
Approximate target

10. Estimating a probability¶

Suppose:

X ~ Normal(0,1)

and we want:

P(X > 1.96)

Simulate:

x = rng.normal(
    size=100000
)

Estimate:

estimate = (
    x > 1.96
).mean()

Because Boolean values act like 0/1 indicators:

True  → 1
False → 0

the mean estimates the probability.

11. Monte Carlo as indicator averaging¶

If:

I(X > 1.96)

is an indicator, then:

E[I(X > 1.96)]
=
P(X > 1.96)

Monte Carlo approximation:

P(X > 1.96)
≈
(1/M) Σ I(X_m > 1.96)

This simple idea is extremely powerful.

12. Estimating expectations¶

Suppose:

X ~ Uniform(0,1)

and we want:

E[X²]

Simulate:

x = rng.uniform(
    0,
    1,
    size=100000
)

Estimate:

np.mean(
    x ** 2
)

Theoretical value:

1/3

Monte Carlo approximates the expectation using a sample average.

13. Monte Carlo integration¶

Many integrals can be written as expectations.

For example:

∫₀¹ x² dx
=
E[X²], where X ~ Uniform(0,1)

Therefore:

x = rng.uniform(
    0,
    1,
    size=M
)

estimate = np.mean(
    x ** 2
)

approximates the integral.

14. Monte Carlo error¶

Suppose:

θ̂_MC

is a Monte Carlo average from M simulated values.

Its standard error typically decreases at approximately:

1 / √M

This means that reducing Monte Carlo error by a factor of 10 generally requires about 100 times as many simulations.

That is an important computational tradeoff.

15. Monte Carlo standard error¶

Suppose we simulate values:

g = x ** 2

Then:

estimate = g.mean()

A Monte Carlo standard error estimate is:

mcse = (
    g.std(ddof=1)
    / np.sqrt(len(g))
)

This quantifies uncertainty from using a finite number of simulations.

16. Simulation size should be justified¶

Avoid choosing:

10,000 simulations

simply because it sounds large.

Ask:

What level of Monte Carlo precision is needed?
How expensive is each simulation?
Is the estimated quantity stable?
What is the Monte Carlo standard error?

Simulation size is part of the analysis design.

17. Visualize convergence¶

Suppose:

draws = rng.normal(
    size=10000
)

Compute cumulative means:

running_mean = (
    np.cumsum(draws)
    / np.arange(
        1,
        len(draws) + 1
    )
)

Plot:

import matplotlib.pyplot as plt

plt.plot(running_mean)
plt.axhline(0)
plt.show()

The running estimate typically stabilizes as more draws accumulate.

18. Vectorized simulation¶

Less efficient:

results = []

for _ in range(10000):
    x = rng.normal()
    results.append(x > 1.96)

Vectorized:

x = rng.normal(
    size=10000
)

estimate = (
    x > 1.96
).mean()

Week 9’s vectorization concepts become especially valuable in simulation.

19. Simulating repeated samples¶

Suppose we want 10,000 samples, each of size 30:

samples = rng.normal(
    loc=0,
    scale=1,
    size=(10000, 30)
)

Then:

sample_means = (
    samples.mean(axis=1)
)

This generates 10,000 sample means with one vectorized operation.

20. Sampling distributions¶

If:

X₁, ..., Xₙ ~ distribution

then repeated samples generate a sampling distribution for a statistic.

Example:

samples = rng.exponential(
    scale=1,
    size=(10000, 30)
)

means = samples.mean(
    axis=1
)

Even though the source distribution is skewed, the distribution of sample means may become approximately normal.

Simulation makes theoretical sampling behavior visible.

21. Inverse-transform sampling¶

Suppose we know a cumulative distribution function:

F(x)

If:

U ~ Uniform(0,1)

then under suitable conditions:

X = F⁻¹(U)

has CDF F.

This is called inverse-transform sampling.

22. Why inverse-transform sampling works¶

Let:

X = F⁻¹(U)

Then:

P(X ≤ x)
=
P(F⁻¹(U) ≤ x)

which implies:

P(U ≤ F(x))
=
F(x)

because:

U ~ Uniform(0,1)

This turns uniform random numbers into samples from another distribution.

23. Exponential example¶

For:

X ~ Exponential(rate = λ)

the CDF is:

F(x) = 1 - exp(-λx)

Set:

u = F(x)

and solve for x:

x = -log(1-u) / λ

Therefore:

u = rng.uniform(
    size=10000
)

lam = 2

x = (
    -np.log(1 - u)
    / lam
)

produces exponential draws.

24. Compare with NumPy¶

NumPy can sample directly:

x_numpy = rng.exponential(
    scale=1 / lam,
    size=10000
)

Compare:

x.mean()
x_numpy.mean()

Both should be close to:

1 / λ

25. Inverse transform and numerical stability¶

Because:

1 - U

is also Uniform(0,1), we can often write:

x = (
    -np.log(u)
    / lam
)

provided:

U ∈ (0,1)

Equivalent formulas may differ numerically at extreme values.

Numerical details matter even in conceptually simple algorithms.

26. Discrete inverse-transform sampling¶

Suppose a discrete variable takes:

A with probability 0.2
B with probability 0.5
C with probability 0.3

Cumulative probabilities:

A: 0.2
B: 0.7
C: 1.0

For:

U ~ Uniform(0,1)

assign:

U < 0.2       → A
0.2 ≤ U < 0.7 → B
0.7 ≤ U       → C

This is the discrete analogue of inverse-transform sampling.

27. Simulating Bernoulli trials¶

A Bernoulli variable can be generated using:

u = rng.uniform(
    size=10000
)

x = (
    u < 0.3
).astype(int)

Then:

x.mean()

should be near:

0.3

NumPy’s direct method:

rng.binomial(
    n=1,
    p=0.3,
    size=10000
)

produces the same distribution.

28. Monte Carlo estimate of π¶

A classic example uses points sampled uniformly in a square.

Generate:

xy = rng.uniform(
    -1,
    1,
    size=(100000, 2)
)

Inside the unit circle:

inside = (
    xy[:, 0] ** 2
    + xy[:, 1] ** 2
    <= 1
)

The area ratio satisfies:

π / 4
≈
proportion inside circle

Therefore:

pi_hat = (
    4 * inside.mean()
)

29. Why the π example is useful¶

Estimating π by simulation is not computationally competitive with direct numerical methods.

Its value is conceptual.

It demonstrates:

random sampling;
geometric probability;
indicator variables;
Monte Carlo averaging;
convergence;
simulation error.

30. Evaluating an estimator by simulation¶

Suppose:

X₁,...,Xₙ ~ Normal(μ, σ²)

and estimator:

μ̂ = sample mean

We can simulate repeated datasets:

samples = rng.normal(
    loc=5,
    scale=2,
    size=(10000, 20)
)

Compute:

estimates = (
    samples.mean(axis=1)
)

Then examine:

estimates.mean()
estimates.std()

Simulation lets us empirically study estimator properties.

31. Bias¶

For estimator:

θ̂

bias is:

Bias(θ̂)
=
E[θ̂] - θ

Simulation estimate:

bias = (
    estimates.mean()
    - true_value
)

Repeated simulation can help assess whether an estimator is approximately unbiased.

32. Variance and MSE¶

Simulation can estimate:

Var(θ̂)

using:

estimates.var(
    ddof=1
)

Mean squared error:

MSE
=
E[(θ̂ - θ)²]

estimated by:

mse = np.mean(
    (
        estimates
        - true_value
    ) ** 2
)

33. Coverage probability¶

Suppose every simulated dataset produces a confidence interval:

[L_m, U_m]

Coverage indicator:

I(L_m ≤ θ ≤ U_m)

Estimated coverage:

coverage = (
    covered.mean()
)

This lets us evaluate whether a nominal:

95%

confidence interval actually covers the true parameter about 95% of the time under the simulation model.

34. Simulation studies are experiments¶

A good simulation study specifies:

data-generating process
parameter values
sample size
number of replications
estimator/procedure
performance metric
random-number strategy

Treat the simulation design as seriously as an empirical experiment.

35. Separate simulation design from execution¶

Instead of hard-coding:

samples = rng.normal(
    5,
    2,
    size=(10000, 20)
)

define:

MU = 5
SIGMA = 2
N = 20
N_SIM = 10000

Then:

samples = rng.normal(
    MU,
    SIGMA,
    size=(N_SIM, N)
)

Named parameters make simulation studies easier to audit.

36. Simulation functions¶

A reusable simulation function might be:

def simulate_mean(
    n,
    mu,
    sigma,
    n_sim,
    rng
):
    samples = rng.normal(
        mu,
        sigma,
        size=(n_sim, n)
    )

    return samples.mean(
        axis=1
    )

Then:

estimates = simulate_mean(
    n=20,
    mu=5,
    sigma=2,
    n_sim=10000,
    rng=rng
)

37. Independent random streams¶

Large simulation projects may need separate random streams.

A simple practical rule is:

Do not repeatedly reset the seed inside a loop.

Bad:

for i in range(100):
    rng = np.random.default_rng(556)
    ...

This regenerates the same sequence repeatedly.

Instead, initialize once and draw sequentially.

38. Parallel simulation requires care¶

If simulations are parallelized, each worker needs an appropriate independent random stream.

Naively giving every worker the same seed can produce duplicated simulations.

This becomes important in later performance-oriented work.

For Week 10, the main principle is:

Random-number management is part of reproducible software design.

39. Validate a simulation¶

Always ask:

Do simulated means match theory?
Are probabilities in [0,1]?
Is the output shape correct?
Are values in the correct support?
Does a histogram resemble the intended distribution?
Does increasing simulation size stabilize the estimate?

Simulation code can run successfully and still simulate the wrong model.

40. Key ideas¶

By the end of Week 10, you should be able to explain:

What simulation means in statistical computing.
What pseudorandom numbers are.
Why random seeds support reproducibility.
Why NumPy generator objects are useful.
How to generate common distributions.
What Monte Carlo estimation is.
How probabilities can be estimated using indicators.
How expectations can be estimated using sample averages.
How Monte Carlo integration works.
Why Monte Carlo error decreases roughly as 1/sqrt(M).
How to estimate Monte Carlo standard error.
Why vectorization is useful for simulation.
What a sampling distribution is.
How inverse-transform sampling works.
How to derive exponential inverse-transform sampling.
How to simulate discrete outcomes from uniforms.
How to evaluate estimators by simulation.
Bias, variance, MSE, and coverage in simulation studies.
Why simulation design parameters should be explicit.
Why random-number management matters in larger computational workflows.

41. Recommended reading¶

NumPy Random Sampling¶

Official NumPy documentation for modern random-number generation.

https://numpy.org/doc/stable/reference/random/index.html

NumPy Generator¶

Reference for the Generator API and supported distributions.

https://numpy.org/doc/stable/reference/random/generator.html

SciPy Statistical Distributions¶

Useful reference for probability distributions, CDFs, inverse CDFs, and random variate generation.

https://docs.scipy.org/doc/scipy/reference/stats.html

Python Data Science Handbook — NumPy¶

Useful background on vectorized numerical operations, which directly support efficient Monte Carlo computation.

https://jakevdp.github.io/PythonDataScienceHandbook/

42. YouTube recommendations¶

1. jbstatistics — “Monte Carlo Simulation”¶

A clear conceptual introduction to Monte Carlo simulation and repeated random sampling. It is useful for reinforcing the idea that simulation approximates probabilities and expectations through repeated experiments.

Recommended use: Watch before the Monte Carlo estimation section.

Find on YouTube

2. Inverse Transform Sampling¶

A focused tutorial on inverse-transform sampling is particularly useful this week because the method connects uniform random-number generation to arbitrary target distributions.

Recommended use: Watch before the inverse-transform sampling exercises.

Find inverse-transform sampling tutorials on YouTube

3. NumPy Random Number Generation¶

A practical NumPy tutorial on default_rng, distributions, seeds, and reproducible random-number generation is useful for the implementation side of the week.

Recommended use: Optional reinforcement before starting the hands-on activity.

Find NumPy RNG tutorials on YouTube

43. Week 10 takeaway¶

The central lesson is:

Simulation converts probability models into computational experiments.

The progression is:

Probability model
      ↓
Random-number generator
      ↓
Simulated samples
      ↓
Statistic / indicator
      ↓
Repeat
      ↓
Monte Carlo approximation
      ↓
Quantify simulation error
      ↓
Statistical insight

Next week we will use repeated simulation ideas to build bootstrap and permutation procedures from scratch.

Why this week matters¶

1. What is simulation?¶

2. Pseudorandom numbers¶

3. Reproducible random-number generation¶

4. Why default_rng()?¶

5. Randomness as a function dependency¶

6. Common distributions¶

7. Verify a simulation¶

8. Simulation variability¶

9. Monte Carlo methods¶

10. Estimating a probability¶

11. Monte Carlo as indicator averaging¶

12. Estimating expectations¶

13. Monte Carlo integration¶

14. Monte Carlo error¶

15. Monte Carlo standard error¶

16. Simulation size should be justified¶

17. Visualize convergence¶

18. Vectorized simulation¶

19. Simulating repeated samples¶

20. Sampling distributions¶

21. Inverse-transform sampling¶

22. Why inverse-transform sampling works¶

23. Exponential example¶

24. Compare with NumPy¶

25. Inverse transform and numerical stability¶

26. Discrete inverse-transform sampling¶

27. Simulating Bernoulli trials¶

28. Monte Carlo estimate of π¶

29. Why the π example is useful¶

30. Evaluating an estimator by simulation¶

31. Bias¶

32. Variance and MSE¶

33. Coverage probability¶

34. Simulation studies are experiments¶

35. Separate simulation design from execution¶

36. Simulation functions¶

37. Independent random streams¶

38. Parallel simulation requires care¶

39. Validate a simulation¶

40. Key ideas¶

41. Recommended reading¶

NumPy Random Sampling¶

NumPy Generator¶

SciPy Statistical Distributions¶

Python Data Science Handbook — NumPy¶

42. YouTube recommendations¶

1. jbstatistics — “Monte Carlo Simulation”¶

2. Inverse Transform Sampling¶

3. NumPy Random Number Generation¶

43. Week 10 takeaway¶

4. Why `default_rng()`?¶