NumPy Random and Simulation

The Generator API

Create one Generator with a seed and draw everything from it.

import numpy as np

rng = np.random.default_rng(seed=42)

print(rng.uniform(0, 1, size=3).round(3))
print(rng.normal(0, 1, size=3).round(3))
print(rng.integers(0, 10, size=5))
print(rng.binomial(n=10, p=0.3, size=5))
print(rng.poisson(lam=3.5, size=5))

$ python generator.py
[0.774 0.439 0.859]
[ 0.941 -1.951 -1.302]
[7 7 7 7 5]
[3 2 5 3 4]
[3 7 1 4 3]

Run it again and you get exactly these numbers — the seed fixes the sequence. integers excludes its upper bound, so 0..9 here; binomial counts successes out of 10 trials, and poisson counts events with an average rate of 3.5.

Sampling, Shuffling and Permuting

import numpy as np

rng = np.random.default_rng(42)
pop = np.array([10, 20, 30, 40, 50])

print(rng.choice(pop, size=3, replace=False))
print(rng.choice(pop, size=5, replace=True))
weighted = rng.choice(pop, size=1000, p=[0.6, 0.1, 0.1, 0.1, 0.1])
print(np.unique(weighted, return_counts=True))

arr = np.arange(8)
perm = rng.permutation(arr)
rng.shuffle(arr)
print(perm, arr)

$ python sampling.py
[50 10 40]
[50 10 40 20 10]
(array([10, 20, 30, 40, 50]), array([605,  93, 106,  98,  98]))
[1 5 4 7 2 6 0 3] [3 5 6 1 7 0 2 4]

Without replacement every pick is distinct; with replacement repeats are allowed. p= weights the picks: over 1,000 draws, 10 (weight 0.6) turns up roughly six times as often as each other value. A handful of draws would be too few to see the weights reliably. permutation returns a shuffled copy, while shuffle reorders arr in place and returns None.

Checking Distributions Against Theory

import numpy as np

rng = np.random.default_rng(42)
N = 100_000
samples = {
    "Normal(5, 2)":      (rng.normal(5, 2, N),          5.0, 2.0),
    "Exponential(2)":    (rng.exponential(2, N),        2.0, 2.0),
    "Poisson(3)":        (rng.poisson(3, N),            3.0, 3 ** 0.5),
    "Binomial(20, 0.3)": (rng.binomial(20, 0.3, N),     6.0, (20 * 0.3 * 0.7) ** 0.5),
}
print(f"{'dist':<18}{'mean':>8}{'theory':>8}{'std':>8}{'theory':>8}")
for name, (x, mu, sd) in samples.items():
    print(f"{name:<18}{x.mean():>8.3f}{mu:>8.3f}{x.std():>8.3f}{sd:>8.3f}")

$ python distributions.py
dist                  mean  theory     std  theory
Normal(5, 2)         4.992   5.000   2.007   2.000
Exponential(2)       1.990   2.000   1.995   2.000
Poisson(3)           2.995   3.000   1.730   1.732
Binomial(20, 0.3)    6.001   6.000   2.044   2.049

With 100,000 draws each sample mean and standard deviation lands within a few thousandths of the textbook value. That comparison is a cheap sanity check whenever you simulate from a distribution you haven’t used before.

Monte Carlo: Estimating π

import numpy as np

rng = np.random.default_rng(42)

for n in [1_000, 100_000, 10_000_000]:
    x = rng.uniform(-1, 1, n)
    y = rng.uniform(-1, 1, n)
    pi_hat = 4 * np.mean(x**2 + y**2 <= 1)
    print(f"n={n:>11,}  pi ~ {pi_hat:.5f}  error {abs(pi_hat - np.pi):.5f}")

$ python monte_carlo_pi.py
n=      1,000  pi ~ 3.08000  error 0.06159
n=    100,000  pi ~ 3.14380  error 0.00221
n= 10,000,000  pi ~ 3.14147  error 0.00013

The fraction of random points inside the unit circle approximates its area, π/4. On average the error shrinks like 1/sqrt(n) — 100 times more points buys about 10 times more accuracy — though any single run scatters around that trend, as these do. That square root is why Monte Carlo is easy to write but expensive to make precise.

Bootstrap Confidence Interval

Resample the data with replacement many times and look at the spread of the statistic.

import numpy as np

rng = np.random.default_rng(42)
ltv = rng.lognormal(mean=5, sigma=1.2, size=500)

idx = rng.integers(0, len(ltv), size=(10_000, len(ltv)))
means = ltv[idx].mean(axis=1)
lo, hi = np.percentile(means, [2.5, 97.5])

print(f"sample mean   {ltv.mean():8.2f}")
print(f"sample median {np.median(ltv):8.2f}")
print(f"95% CI mean   [{lo:.2f}, {hi:.2f}]")

$ python bootstrap.py
sample mean     288.81
sample median   148.99
95% CI mean   [249.79, 332.97]

idx holds 10,000 rows of 500 random positions, so ltv[idx] is 10,000 resamples at once and .mean(axis=1) reduces each row. The middle 95% of those means is the interval. The mean sits well above the median because log-normal revenue data has a long right tail.

Loop vs apply_along_axis vs a Real Reduction

import timeit
import numpy as np

rng = np.random.default_rng(42)
data = rng.lognormal(5, 1.2, 500)
idx = rng.integers(0, 500, size=(5_000, 500))

loop = lambda: np.array([data[row].mean() for row in idx])
applied = lambda: np.apply_along_axis(lambda r: data[r].mean(), 1, idx)
reduced = lambda: data[idx].mean(axis=1)

for name, f in [("python loop", loop), ("apply_along_axis", applied),
                ("mean(axis=1)", reduced)]:
    print(f"{name:17s} {min(timeit.repeat(f, number=1, repeat=5)):.4f} s")
print(np.allclose(loop(), applied()), np.allclose(loop(), reduced()))

$ python bootstrap_speed.py
python loop       0.0155 s
apply_along_axis  0.0158 s
mean(axis=1)      0.0056 s
True True

All three compute the same 5,000 means. apply_along_axis still calls a Python function once per row, so taking the best of five runs it matches the plain loop almost exactly. Only mean(axis=1) moves the work into one compiled reduction, and it is roughly three times faster here. Absolute timings vary between machines.

Independent Streams with SeedSequence

import numpy as np

children = np.random.SeedSequence(2024).spawn(3)
rngs = [np.random.default_rng(c) for c in children]
for i, r in enumerate(rngs):
    print(i, r.integers(0, 100, size=4))

again = [np.random.default_rng(c) for c in np.random.SeedSequence(2024).spawn(3)]
print(again[1].integers(0, 100, size=4))

$ python spawn.py
0 [67 65 15 52]
1 [27 12 37 45]
2 [76 92 82 12]
[27 12 37 45]

Each child seed gives its own stream, suitable for one worker or one experiment each, and the whole family is reproducible from the single base seed — the last line repeats stream 1 exactly. Deriving seeds by hand (seed + i * 1000) gives no such independence guarantee.

Practice

1. Simulate 100,000 rolls of two dice and print the probability of rolling a total of 7 (theory: 1/6).
2. Bootstrap a 90% confidence interval for the median of the ltv data above.
3. Show that rng.shuffle returns None.

import numpy as np

rng = np.random.default_rng(7)
rolls = rng.integers(1, 7, size=(100_000, 2)).sum(axis=1)
print(f"P(7) ~ {np.mean(rolls == 7):.4f}  theory {1/6:.4f}")

$ python practice_dice.py
P(7) ~ 0.1666  theory 0.1667

import numpy as np

rng = np.random.default_rng(42)
ltv = rng.lognormal(mean=5, sigma=1.2, size=500)
idx = rng.integers(0, 500, size=(10_000, 500))
medians = np.median(ltv[idx], axis=1)
print(np.percentile(medians, [5, 95]).round(2))

$ python practice_median_ci.py
[129.17 165.76]

import numpy as np

rng = np.random.default_rng(0)
arr = np.arange(5)
print(rng.shuffle(arr))
print(arr)

$ python practice_shuffle.py
None
[2 4 3 0 1]