Activity 5.2 - STA556 Statistics and Data Science Computing Workflows

This page is an ipython notebook, meaning you can download it and run the executable cells. It directly follows Activity 5.2. Reference answers appear after the questions, challenges, predictions, and reflections. They show one reasonable solution or interpretation, but other correct approaches may exist.

Random Number Generation, Inverse Transform Sampling & Monte Carlo¶

Tools: Python, NumPy, matplotlib, Jupyter/VS Code

Learning objectives¶

By the end of this tutorial, you should be able to:

use NumPy’s modern random-number generator;
reproduce simulations using a fixed seed;
simulate from common distributions;
verify simulation output against theory;
estimate probabilities and expectations with Monte Carlo;
calculate Monte Carlo standard errors;
visualize convergence;
vectorize repeated simulation;
generate sampling distributions;
implement inverse-transform sampling;
simulate discrete variables from uniforms;
evaluate estimators using repeated simulation;
estimate bias, variance, MSE, and coverage.

Part 0 — Set up¶

Create:

notebooks/week10_simulation.ipynb

Import:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Create:

rng = np.random.default_rng(556)

Part 1 — Reproducible RNG¶

Run:

rng = np.random.default_rng(556)

rng.normal(
    size=5
)

Restart:

rng = np.random.default_rng(556)

rng.normal(
    size=5
)

Question¶

Are the values identical?

Why?

Part 2 — Different seeds¶

Compare:

rng1 = np.random.default_rng(556)
rng2 = np.random.default_rng(557)

Generate:

rng1.normal(size=5)
rng2.normal(size=5)

Reflection¶

What role does the seed play?

Reference answer¶

Recreating a generator with the same seed and requesting the same sequence of draws reproduces the same values. Different seeds produce different pseudorandom sequences.

Part 3 — Simulate common distributions¶

Generate 10,000 draws from:

Normal¶

normal = rng.normal(
    10,
    2,
    size=10000
)

Uniform¶

uniform = rng.uniform(
    0,
    1,
    size=10000
)

Binomial¶

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

Poisson¶

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

Reference answer¶

The seed initializes the pseudorandom generator, making the sequence reproducible.

Part 4 — Verify theoretical moments¶

For the normal sample:

normal.mean()
normal.std()

For the uniform sample:

uniform.mean()
uniform.var()

For the Poisson sample:

poisson.mean()
poisson.var()

Questions¶

How close are sample moments to theoretical values?
Why are they not exactly equal?

Part 5 — Visualize distributions¶

Create histograms for the four simulated variables.

Example:

plt.hist(
    normal,
    bins=30
)

plt.show()

Reflection¶

Does the visual shape match the intended model?

Reference answer¶

Sample moments should be close to their theoretical values but not exactly equal because finite samples contain random variation.

Part 6 — Randomness as a dependency¶

Write:

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

Call:

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

Question¶

Why is passing rng preferable to creating a new seed inside the function?

Reference answer¶

The histograms should resemble the intended theoretical distributions in broad shape, but random sampling prevents a perfect match.

Part 7 — Monte Carlo probability¶

Estimate:

P(Z > 1.96)

where:

Z ~ Normal(0,1)

Use:

z = rng.normal(
    size=100000
)

estimate = (
    z > 1.96
).mean()

estimate

Compare with the known probability, approximately 0.025.

Reference answer¶

Passing rng makes randomness an explicit dependency and avoids silently reseeding inside the function.

Part 8 — Indicator interpretation¶

Inspect:

(z > 1.96)[:10]

Then:

(z > 1.96).astype(int)[:10]

Explain why the mean of the Boolean values estimates a probability.

Part 9 — Estimate an expectation¶

Estimate:

E[X²]

for:

X ~ Uniform(0,1)

Use:

x = rng.uniform(
    size=100000
)

estimate = np.mean(
    x ** 2
)

Compare with:

1/3

Reference answer¶

Boolean True/False values behave like 1/0 when averaged, so their mean is the proportion of simulated draws for which the event occurred.

Part 10 — Monte Carlo standard error¶

Set:

g = x ** 2

Compute:

estimate = g.mean()

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

estimate
mcse

Question¶

What uncertainty does mcse represent?

Part 11 — Effect of simulation size¶

Repeat the same expectation estimate using:

M = 100
M = 1,000
M = 10,000
M = 100,000

For each:

estimate the expectation;
calculate MCSE.

Store results in a DataFrame.

Reflection¶

How does MCSE change as M increases?

Reference answer¶

MCSE measures error from using a finite number of Monte Carlo replications. For a probability it is approximately sqrt(p_hat*(1-p_hat)/M).

Part 12 — Running convergence¶

Generate:

draws = rng.normal(
    size=10000
)

Compute:

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

Plot:

plt.plot(running_mean)
plt.axhline(0)
plt.xlabel("Number of draws")
plt.ylabel("Running mean")
plt.show()

Question¶

What does the plot show?

Reference answer¶

MCSE decreases at about 1/sqrt(M), so halving MCSE requires roughly four times as many simulations.

Part 13 — Loop vs. vectorized simulation¶

Loop:

results = []

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

Vectorized:

values = rng.normal(
    size=10000
)

vectorized_result = (
    values > 1.96
).mean()

Compare the implementations.

Reference answer¶

The running estimate should fluctuate substantially early and then stabilize around the target as more simulations accumulate.

Part 14 — Benchmark¶

Use:

%timeit

to compare loop and vectorized versions.

Reflection¶

How does Week 9’s vectorization material improve simulation workflows?

Reference answer¶

The vectorized implementation is generally faster because it moves repeated work into optimized NumPy routines.

Part 15 — Repeated samples¶

Generate:

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

Inspect:

samples.shape

Calculate:

means = samples.mean(
    axis=1
)

Inspect:

means.shape

Question¶

What does each row represent?

Reference answer¶

Vectorization makes repeated simulation and summaries across many samples both faster and often clearer.

Part 16 — Sampling distribution of the mean¶

Plot:

plt.hist(
    means,
    bins=30
)

plt.show()

Check:

means.mean()
means.std()

Compare with the theoretical standard error:

1 / sqrt(30)

Reference answer¶

Each row represents one independently simulated sample. The row means form a sampling distribution whose standard deviation should be close to the theoretical standard error.

Part 17 — CLT experiment¶

Generate exponential samples:

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

means = samples.mean(
    axis=1
)

Plot:

plt.hist(
    means,
    bins=30
)

plt.show()

Reflection¶

How does the sampling distribution compare with the source distribution?

Part 18 — Inverse transform: exponential¶

Generate uniforms:

u = rng.uniform(
    size=100000
)

Set:

lam = 2

Transform:

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

Check:

x_inverse.mean()

Compare with:

1 / lam

Reference answer¶

The sampling distribution of the mean is more concentrated than the source distribution and becomes more nearly normal as sample size grows.

Part 19 — Compare with direct exponential sampling¶

Generate:

x_direct = rng.exponential(
    scale=1 / lam,
    size=100000
)

Compare:

x_inverse.mean()
x_direct.mean()

and:

x_inverse.var()
x_direct.var()

Plot both distributions.

Reference answer¶

Inverse-transform and direct exponential samples should have similar distributions/moments but not identical individual draws.

Part 20 — Verify the inverse transform¶

For selected probabilities:

probs = np.array([
    0.1,
    0.5,
    0.9
])

Compute:

quantiles = (
    -np.log(1 - probs)
    / lam
)

Interpret the result as theoretical exponential quantiles.

Part 21 — Bernoulli from uniforms¶

Generate:

u = rng.uniform(
    size=10000
)

Create:

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

Check:

bernoulli.mean()

Compare with:

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

direct.mean()

Reference answer¶

The Bernoulli sample mean estimates the success probability. Inverse-CDF transformation maps uniform draws into draws from the target distribution.

Part 22 — Discrete inverse-transform sampling¶

Suppose:

P(A)=0.2
P(B)=0.5
P(C)=0.3

Generate:

u = rng.uniform(
    size=10000
)

Assign:

outcome = np.where(
    u < 0.2,
    "A",
    np.where(
        u < 0.7,
        "B",
        "C"
    )
)

Check proportions:

pd.Series(
    outcome
).value_counts(
    normalize=True
)

Part 23 — Estimate π¶

Generate:

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

Determine:

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

Estimate:

pi_hat = (
    4 * inside.mean()
)

pi_hat

Compare with:

np.pi

Reference answer¶

The simulation estimates π using the probability that a uniform point in a square lies inside an inscribed circle; multiply the estimated fraction by 4.

Part 24 — Visualize the π experiment¶

Use a smaller sample:

xy_small = rng.uniform(
    -1,
    1,
    size=(2000, 2)
)

Create an indicator for points inside the circle.

Plot the points with separate markers or categories for inside/outside.

Question¶

What geometric probability is being estimated?

Part 25 — Monte Carlo integration challenge¶

Estimate:

∫₀¹ exp(-x²) dx

using:

X ~ Uniform(0,1)

Your estimate is:

x = rng.uniform(
    size=100000
)

estimate = np.mean(
    np.exp(
        -(x ** 2)
    )
)

Calculate an MCSE.

Reference answer¶

Monte Carlo integration estimates an integral using an average of function values over random draws plus the appropriate scaling factor.

Part 26 — Simulate an estimator¶

Set:

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

Generate:

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

Estimate:

mu_hat = samples.mean(
    axis=1
)

Inspect:

mu_hat.mean()
mu_hat.std()

Part 27 — Bias¶

Calculate:

bias = (
    mu_hat.mean()
    - MU
)

bias

Question¶

Does the sample mean appear unbiased?

Part 28 — Variance and MSE¶

Compute:

variance = mu_hat.var(
    ddof=1
)

and:

mse = np.mean(
    (
        mu_hat - MU
    ) ** 2
)

Compare:

variance
mse

Why are they close when bias is close to zero?

Reference answer¶

Bias measures systematic error, variance measures estimator variability, and MSE combines them: approximately variance + bias**2.

Part 29 — Compare two estimators¶

For each simulated sample calculate:

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

and:

median_estimator = (
    np.median(
        samples,
        axis=1
    )
)

Compare:

bias;
variance;
MSE.

Reflection¶

Which estimator performs better under normal data for estimating the population mean?

Reference answer¶

If bias is near zero, MSE is close to variance. Under normal data, the sample mean is typically more efficient than the sample median for estimating the population mean.

Part 30 — Confidence interval coverage¶

For each sample compute:

means = samples.mean(
    axis=1
)

ses = samples.std(
    axis=1,
    ddof=1
) / np.sqrt(N)

Create approximate 95% intervals:

lower = means - 1.96 * ses
upper = means + 1.96 * ses

Check:

covered = (
    (lower <= MU)
    & (MU <= upper)
)

Estimate:

coverage = (
    covered.mean()
)

coverage

Reference answer¶

Coverage is the fraction of simulated intervals containing the true parameter. A nominal 95% method should have coverage near 0.95 over many replications.

Part 31 — Simulation function¶

Write:

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

    return samples.mean(
        axis=1
    )

Test with:

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

Reference answer¶

A reusable simulation function should accept scenario parameters and an RNG explicitly, and return simulated results rather than only printing them.

Part 32 — Compare sample sizes¶

Run the function for:

n = 5
n = 20
n = 100

For each calculate:

bias;
standard deviation;
MSE.

Create a summary DataFrame.

Question¶

How does estimator variability change with sample size?

Reference answer¶

Estimator variability decreases as sample size increases, commonly at a rate related to 1/sqrt(n).

Part 33 — Avoid resetting seeds inside loops¶

Compare this bad pattern:

results = []

for _ in range(5):
    rng_bad = np.random.default_rng(556)

    results.append(
        rng_bad.normal()
    )

Inspect:

results

Now initialize once:

rng_good = np.random.default_rng(556)

results = [
    rng_good.normal()
    for _ in range(5)
]

Question¶

Why is the first pattern wrong?

Reference answer¶

Resetting the same seed inside a loop can generate repeated random streams and destroy the intended independence of replications.

Part 34 — Simulation study design¶

Write a study plan before coding.

For example:

Goal:
Data-generating distribution:
True parameter:
Sample size:
Number of simulations:
Estimator:
Performance measures:
Random seed:

Use your own design.

Then implement it.

Reference answer¶

A simulation plan should state the data-generating process, parameter values, sample sizes, number of replications, estimators, metrics, RNG strategy, and validation checks.

Part 35 — Mini simulation challenge¶

Investigate estimation of a population mean under a skewed distribution.

Use:

X ~ Exponential(mean=1)

Compare:

sample mean
sample median

for estimating:

population mean = 1

Use:

n = 10, 30, 100

and:

10,000 simulations

For each estimator and sample size, report:

bias;
variance;
MSE.

Reflection¶

Does the median estimate the same population parameter as the mean?

This is an important conceptual check.

Reference answer¶

Compare estimators using bias, variance, and MSE across scenarios, and interpret differences in terms of sample size and the generating distribution.

Part 36 — Validate the simulation¶

For your mini study, verify:

samples.shape

Check:

np.isfinite(
    samples
).all()

Compare simulated moments with theoretical moments.

Check that probabilities are between:

0 and 1

Principle¶

A simulation that runs is not necessarily a simulation of the intended model.

Reference answer¶

Validate against known theoretical moments/probabilities where possible, check finite outputs and valid ranges, and confirm the RNG design behaves as intended.

Part 37 — Git checkpoint¶

Run:

git status

Stage:

git add .

Commit:

git commit -m "Complete Week 10 simulation and Monte Carlo exercises"

Push:

git push

Part 38 — Final reflection¶

Answer in Markdown.

1. Randomness¶

What does a random seed control?

2. Generator design¶

Why pass an RNG object into a simulation function?

3. Monte Carlo¶

What is the central idea of Monte Carlo estimation?

4. Probability¶

Why does the mean of an indicator estimate a probability?

5. MCSE¶

What does Monte Carlo standard error measure?

6. Vectorization¶

Why is vectorization especially useful for simulation?

7. Sampling distributions¶

What does a repeated-sampling simulation approximate?

8. Inverse transform¶

How can Uniform(0,1) draws be transformed into another distribution?

9. Estimator evaluation¶

How do bias, variance, and MSE differ?

10. Reproducibility¶

Why is random-number management part of scientific reproducibility?

Completion checklist¶

Created Week 10 notebook
Used np.random.default_rng()
Reproduced results with a fixed seed
Compared different seeds
Simulated Normal, Uniform, Binomial, and Poisson variables
Checked simulated moments against theory
Visualized simulated distributions
Passed an RNG object into a function
Estimated a probability with Monte Carlo
Interpreted Boolean indicators
Estimated an expectation
Calculated Monte Carlo standard error
Compared simulation sizes
Visualized running convergence
Compared loop and vectorized simulation
Simulated repeated samples in a 2D array
Generated a sampling distribution
Explored the CLT with skewed data
Implemented exponential inverse-transform sampling
Compared inverse-transform and direct sampling
Simulated Bernoulli outcomes from uniforms
Implemented a discrete inverse transform
Estimated π by Monte Carlo
Performed Monte Carlo integration
Simulated an estimator
Calculated bias
Calculated variance
Calculated MSE
Compared two estimators
Estimated confidence-interval coverage
Wrote a reusable simulation function
Compared estimator behavior across sample sizes
Demonstrated why seeds should not be reset inside loops
Designed a simulation study before coding
Completed the mini simulation challenge
Validated simulation output
Committed Week 10 work to Git
Pushed work to GitHub

What you should now understand¶

Probability model
      ↓
Pseudorandom generator
      ↓
Simulated samples
      ↓
Vectorized statistics
      ↓
Repeated experiments
      ↓
Monte Carlo estimate
      ↓
Simulation uncertainty
      ↓
Evaluate statistical procedures

Next week we will build bootstrap and permutation procedures from scratch, using the repeated-sampling framework developed here.

Reference answer¶

Seeds control reproducibility; RNG objects make randomness explicit. Monte Carlo uses simulated averages to estimate expectations/probabilities, MCSE quantifies simulation error, and validation is essential.