NumPy Random and Simulation

Real-World Scenario

A data scientist runs a Monte Carlo simulation to estimate the 95% confidence interval for customer lifetime value. Without a fixed seed, they can’t reproduce the result. With np.random.default_rng(42), the exact simulation is reproducible, and bootstrapping gives honest confidence intervals without assumptions about the underlying distribution.

The Generator API (Modern NumPy)

import numpy as np

# Create a reproducible Generator
rng = np.random.default_rng(seed=42)

# Basic distributions
uniform  = rng.uniform(low=0, high=1, size=5)
normal   = rng.normal(loc=0, scale=1, size=5)
integers = rng.integers(low=0, high=10, size=5)   # [low, high)
binomial = rng.binomial(n=10, p=0.3, size=5)      # 10 trials, 30% success
poisson  = rng.poisson(lam=3.5, size=5)           # average 3.5 events

print(f"uniform:  {uniform.round(3)}")
print(f"normal:   {normal.round(3)}")
print(f"integers: {integers}")
print(f"binomial: {binomial}")
print(f"poisson:  {poisson}")

# Sampling from arrays
population = np.array([10, 20, 30, 40, 50])

without_replace = rng.choice(population, size=3, replace=False)
with_replace    = rng.choice(population, size=5, replace=True)
weighted        = rng.choice(population, size=3, p=[0.5, 0.2, 0.1, 0.1, 0.1])
print(f"\nwithout replacement: {without_replace}")
print(f"with replacement:    {with_replace}")
print(f"weighted:            {weighted}")

# Shuffle in-place
arr = np.arange(10)
rng.shuffle(arr)
print(f"\nshuffled: {arr}")

# Permutation (returns new array, original unchanged)
original  = np.arange(10)
permuted  = rng.permutation(original)
print(f"original:  {original}")
print(f"permuted:  {permuted}")

Common Distributions

import numpy as np

rng = np.random.default_rng(42)
N   = 100_000

distributions = {
    # Continuous
    "Uniform(0,1)":           rng.uniform(0, 1, N),
    "Normal(0,1)":            rng.standard_normal(N),
    "Normal(mean=5, std=2)":  rng.normal(5, 2, N),
    "Exponential(scale=2)":   rng.exponential(scale=2, size=N),
    "Log-Normal":             rng.lognormal(mean=0, sigma=1, size=N),
    "Beta(2,5)":              rng.beta(2, 5, N),
    "Gamma(shape=2,scale=2)": rng.gamma(shape=2, scale=2, size=N),

    # Discrete
    "Poisson(lam=3)":         rng.poisson(lam=3, size=N),
    "Binomial(n=20, p=0.3)":  rng.binomial(n=20, p=0.3, size=N),
    "Geometric(p=0.2)":       rng.geometric(p=0.2, size=N),
}

print(f"{'Distribution':<30} {'Mean':>8} {'Std':>8} {'Min':>8} {'Max':>8}")
print("-" * 66)
for name, samples in distributions.items():
    print(f"{name:<30} {samples.mean():>8.3f} {samples.std():>8.3f} "
          f"{samples.min():>8.3f} {samples.max():>8.3f}")

Monte Carlo Simulation

import numpy as np

rng = np.random.default_rng(42)

# ── Example 1: Estimate π using Monte Carlo ────────────────────────
def estimate_pi(n_samples: int, rng: np.random.Generator) -> float:
    """Drop n_samples points in a unit square, count those inside the unit circle."""
    x, y  = rng.uniform(-1, 1, n_samples), rng.uniform(-1, 1, n_samples)
    inside = (x**2 + y**2) <= 1
    return 4 * inside.mean()

for n in [1_000, 10_000, 100_000, 1_000_000]:
    pi_est = estimate_pi(n, rng)
    print(f"n={n:>10,}  π ≈ {pi_est:.6f}  error = {abs(pi_est - np.pi):.6f}")

# ── Example 2: Portfolio Value at Risk (VaR) ──────────────────────
# Simulate daily returns for a portfolio of 3 assets
N_DAYS      = 252      # trading days per year
N_SIMS      = 50_000
INITIAL_VAL = 100_000  # $100k portfolio

# Asset parameters: [tech_stock, bond_etf, gold_etf]
daily_returns_mean = np.array([0.0008, 0.0002, 0.0003])
daily_returns_std  = np.array([0.015,  0.005,  0.010])
weights            = np.array([0.5, 0.3, 0.2])   # portfolio weights

# Simulate N_SIMS paths × N_DAYS returns
asset_returns = rng.normal(
    loc=daily_returns_mean,
    scale=daily_returns_std,
    size=(N_SIMS, N_DAYS, 3),
)
portfolio_daily = (asset_returns * weights).sum(axis=2)  # (N_SIMS, N_DAYS)
portfolio_path  = (1 + portfolio_daily).cumprod(axis=1)  # (N_SIMS, N_DAYS)
final_values    = INITIAL_VAL * portfolio_path[:, -1]    # (N_SIMS,)

# Value at Risk and Expected Shortfall
VaR_95  = np.percentile(final_values, 5)
VaR_99  = np.percentile(final_values, 1)
ES_95   = final_values[final_values <= VaR_95].mean()   # average loss beyond VaR

print(f"\nPortfolio Monte Carlo ({N_SIMS:,} simulations, {N_DAYS} days):")
print(f"  Initial value:         ${INITIAL_VAL:>12,.0f}")
print(f"  Expected final value:  ${final_values.mean():>12,.0f}")
print(f"  95% VaR:               ${VaR_95:>12,.0f}  (5% chance of losing more)")
print(f"  99% VaR:               ${VaR_99:>12,.0f}  (1% chance of losing more)")
print(f"  Expected Shortfall:    ${ES_95:>12,.0f}  (average loss if < 5% VaR)")

Bootstrap Confidence Intervals

import numpy as np

rng = np.random.default_rng(42)

def bootstrap_ci(
    data:       np.ndarray,
    statistic:  callable,
    n_resamples: int = 10_000,
    ci:         float = 0.95,
) -> tuple[float, float, float]:
    """
    Compute a bootstrap confidence interval for any statistic.
    Returns (point_estimate, lower_bound, upper_bound).
    """
    n     = len(data)
    stats = np.empty(n_resamples)

    for i in range(n_resamples):
        resample  = rng.choice(data, size=n, replace=True)
        stats[i]  = statistic(resample)

    alpha = 1 - ci
    lower = np.percentile(stats, 100 * alpha / 2)
    upper = np.percentile(stats, 100 * (1 - alpha / 2))
    return statistic(data), lower, upper


# Simulate customer LTV data (log-normal, as typical revenue data)
ltv_data = rng.lognormal(mean=5, sigma=1.2, size=500)

# Bootstrap CIs for mean, median, and 90th percentile
for stat_name, stat_fn in [
    ("Mean",       np.mean),
    ("Median",     np.median),
    ("90th pctile", lambda x: np.percentile(x, 90)),
]:
    est, lo, hi = bootstrap_ci(ltv_data, stat_fn)
    print(f"{stat_name:<15}: ${est:>8.2f}  95% CI: [${lo:>8.2f}, ${hi:>8.2f}]")

# Vectorized bootstrap (much faster for large n_resamples)
def bootstrap_ci_vectorized(
    data:       np.ndarray,
    statistic:  callable,
    n_resamples: int = 10_000,
    ci:         float = 0.95,
) -> tuple[float, float, float]:
    """Vectorized version — 10-100x faster than the loop version."""
    n        = len(data)
    indices  = rng.integers(0, n, size=(n_resamples, n))
    resamples = data[indices]                    # (n_resamples, n)
    stats    = np.apply_along_axis(statistic, 1, resamples)
    alpha    = 1 - ci
    return (
        statistic(data),
        np.percentile(stats, 100 * alpha / 2),
        np.percentile(stats, 100 * (1 - alpha / 2)),
    )

import timeit
t_loop = timeit.timeit(lambda: bootstrap_ci(ltv_data, np.mean, 5000), number=1)
t_vec  = timeit.timeit(lambda: bootstrap_ci_vectorized(ltv_data, np.mean, 5000), number=1)
print(f"\nBootstrap speed (5000 resamples):")
print(f"  Loop:       {t_loop:.3f}s")
print(f"  Vectorized: {t_vec:.3f}s  ({t_loop/t_vec:.0f}x speedup)")

Reproducible Experiment Pattern

import numpy as np

def make_experiment_rng(base_seed: int, experiment_id: int) -> np.random.Generator:
    """
    Create a deterministic but independent Generator for each experiment.
    Avoids correlated random streams between experiments.
    """
    return np.random.default_rng(base_seed + experiment_id * 1000)


def run_experiment(seed: int, n_samples: int = 1000) -> dict:
    """A fully reproducible experiment."""
    rng = np.random.default_rng(seed)

    # Data generation
    X = rng.normal(0, 1, (n_samples, 5))
    y = (X[:, 0] * 2 + X[:, 1] - X[:, 2] + rng.normal(0, 0.5, n_samples)) > 0

    # Train/test split (reproducible)
    indices    = rng.permutation(n_samples)
    split      = int(n_samples * 0.8)
    train_idx  = indices[:split]
    test_idx   = indices[split:]

    return {
        "seed":       seed,
        "n_train":    len(train_idx),
        "n_test":     len(test_idx),
        "pos_rate":   round(float(y[train_idx].mean()), 4),
        "X_mean":     round(float(X[train_idx].mean()), 4),
    }


# Running the same experiment twice gives identical results
result_a = run_experiment(seed=42)
result_b = run_experiment(seed=42)
assert result_a == result_b
print("Same seed → identical results ✓")

# Different seeds give different results
result_c = run_experiment(seed=99)
print(f"Seed 42 pos_rate: {result_a['pos_rate']}")
print(f"Seed 99 pos_rate: {result_c['pos_rate']}")

# Parallel experiments with independent streams
results = [run_experiment(seed=i) for i in range(5)]
pos_rates = [r["pos_rate"] for r in results]
print(f"\nPos rates across 5 seeds: {pos_rates}")
print(f"Variation: std={np.std(pos_rates):.4f}")