64. Sampling

Most of this course has assumed you can compute with distributions analytically. Beyond the textbook cases you usually can't — and the practical answer is to draw samples and compute with those instead.

Everything in this section rests on one starting point: a stream of independent Uniform(0,1) numbers. This lesson turns that stream into anything else.

Uniform random numbers

Computers produce pseudo-random numbers: deterministic sequences that pass statistical tests for randomness. Modern generators like the Mersenne Twister and PCG64 have enormous periods and good equidistribution.

Determinism is a feature, not a flaw. Seeding makes results reproducible, which is why every simulation in this course sets one.

Pseudo-random generators are not cryptographically secure — observing enough output reveals the internal state. Use secrets or os.urandom when unpredictability matters, and a seeded PRNG when reproducibility does.

Inverse transform sampling

The general method, and the one you met twice already.

Claim. If U \sim \text{Uniform}(0,1) and F is a CDF, then X = F^{-1}(U) has CDF F.

Proof, using that F is non-decreasing:

P(X \le x) = P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x)

That's it — one line, and it works for any distribution whose CDF you can invert.

This is the probability integral transform from §2, run backwards. Forwards it turns any continuous variable into a uniform; backwards it turns a uniform into anything.

Example: the Exponential. F(x) = 1 - e^{-\lambda x}, so solving u = 1 - e^{-\lambda x} gives

X = -\frac{\ln(1-U)}{\lambda} = -\frac{\ln U}{\lambda}

(the last step because 1 - U is also uniform). We derived this in §2 as a Jacobian exercise; here it's an instance of a general recipe.

Discrete version. Cumulate the PMF and find the first level exceeding U — which is what np.searchsorted on a cumulative sum does.

The limitation: you need F^{-1} in closed form. The Normal doesn't have one, which is why it needs its own method.

Sampling the Normal

Box–Muller. From two independent uniforms:

Z_1 = \sqrt{-2\ln U_1}\cos(2\pi U_2), \qquad Z_2 = \sqrt{-2\ln U_1}\sin(2\pi U_2)

gives two independent standard Normals. The trick is polar coordinates: the squared radius of a 2-D standard Normal is Exponential (which is invertible), and the angle is uniform.

Ziggurat. What libraries actually use — a rejection method with precomputed layers, several times faster.

Once you have a standard Normal, any other follows: X = \mu + \sigma Z.

Other useful constructions

Many distributions are cheapest to sample by building them from their definition rather than inverting anything:

Target	Construction
Binomial(n,p)	sum of n Bernoullis
Gamma(k,\lambda), integer k	sum of k Exponentials
\chi^2_k	sum of k squared standard Normals
Beta(a,b)	\frac{G_1}{G_1+G_2} with G_i \sim Gamma
Student's t_k	Z/\sqrt{\chi^2_k/k}
Poisson(\lambda)	count Exponential gaps until they exceed 1

These are exactly the relationships from §1 and §2, now read as algorithms.

Sampling without a normalising constant

The hardest and most important case: you know p(x) \propto \tilde p(x) but cannot compute the normalising constant \int \tilde p.

This is the standard situation in Bayesian inference. The posterior is p(\theta \mid \text{data}) \propto p(\text{data} \mid \theta)p(\theta), and the denominator is an integral over the whole parameter space — usually intractable (§7).

Three responses, in increasing generality:

• Rejection sampling — next lesson.
• Importance sampling — the last lesson of this course.
• MCMC — build a Markov chain whose stationary distribution is p (§8's detailed balance), then run it. This is what makes modern Bayesian computation possible.

Worked example

Sample from a distribution with density f(x) = 3x^2 on [0,1].

CDF: F(x) = x^3 (computed in §1).

Invert: u = x^3 \implies x = u^{1/3}.

Algorithm: draw U \sim \text{Uniform}(0,1), return U^{1/3}.

Check it. The mean should be

E[X] = \int_0^1 x \cdot 3x^2\,dx = \frac{3}{4} = 0.75

and the sampler gives 0.75. The proportion below 0.5 should be F(0.5) = 0.125, and it is.

Notice how little work this took: one integration to get F, one inversion, and the sampler is a single expression. When F^{-1} exists in closed form, inverse transform is unbeatable.

Doing it in Python

Inverse transform for the Exponential, verified:

import numpy as np
from scipy.stats import kstest, expon

rng = np.random.default_rng(0)
lam, n = 2.0, 200_000

U = rng.uniform(0, 1, n)
X = -np.log(U) / lam                 # inverse transform

print(f"mean     : {X.mean():.4f}   theory {1/lam:.4f}")
print(f"variance : {X.var():.4f}   theory {1/lam**2:.4f}")
print(f"KS test against Exponential(2): p = {kstest(X, expon(scale=1/lam).cdf).pvalue:.4f}")
print("\n(A large p-value means we cannot distinguish our samples from the real thing.)")

The 3x^2 example, and the general recipe:

import numpy as np

rng = np.random.default_rng(1)
n = 500_000

U = rng.uniform(0, 1, n)
X = U ** (1/3)                       # F(x) = x^3, so F^-1(u) = u^(1/3)

print(f"mean          : {X.mean():.4f}   theory 0.7500")
print(f"P(X < 0.5)    : {(X < 0.5).mean():.4f}   theory {0.5**3:.4f}")
print(f"P(X < 0.9)    : {(X < 0.9).mean():.4f}   theory {0.9**3:.4f}")
print(f"median        : {np.median(X):.4f}   theory {0.5**(1/3):.4f}")

Discrete inverse transform — cumulate and search:

import numpy as np

rng = np.random.default_rng(2)

values = np.array([1, 2, 3, 4, 5])
probs = np.array([0.1, 0.2, 0.4, 0.2, 0.1])
cdf = np.cumsum(probs)

n = 400_000
U = rng.uniform(0, 1, n)
X = values[np.searchsorted(cdf, U)]

print(f"{'value':>7} {'target':>9} {'sampled':>9}")
for v, p in zip(values, probs):
    print(f"{v:>7} {p:>9.4f} {(X == v).mean():>9.4f}")

Box–Muller from scratch:

import numpy as np
from scipy.stats import kstest, norm

rng = np.random.default_rng(3)
n = 200_000

U1, U2 = rng.uniform(0, 1, n), rng.uniform(0, 1, n)
R = np.sqrt(-2 * np.log(U1))
Z1 = R * np.cos(2 * np.pi * U2)
Z2 = R * np.sin(2 * np.pi * U2)

print(f"Z1: mean {Z1.mean():+.4f}  sd {Z1.std():.4f}")
print(f"Z2: mean {Z2.mean():+.4f}  sd {Z2.std():.4f}")
print(f"correlation between Z1 and Z2: {np.corrcoef(Z1, Z2)[0,1]:+.4f}  (independent)")
print(f"\nKS test on Z1 against N(0,1): p = {kstest(Z1, 'norm').pvalue:.4f}")

Building distributions from their definitions:

import numpy as np
from scipy.stats import chi2, beta as beta_dist, t as t_dist

rng = np.random.default_rng(4)
n = 200_000

# chi-squared: sum of k squared standard normals
k = 5
X = (rng.standard_normal((n, k)) ** 2).sum(axis=1)
print(f"chi2({k}) : mean {X.mean():.4f} (theory {chi2.mean(k)})  "
      f"var {X.var():.4f} (theory {chi2.var(k)})")

# Beta: G1 / (G1 + G2) with gammas
a, b = 2.0, 5.0
G1, G2 = rng.gamma(a, 1, n), rng.gamma(b, 1, n)
B = G1 / (G1 + G2)
print(f"Beta({a},{b}): mean {B.mean():.4f} (theory {beta_dist.mean(a, b):.4f})")

# Student's t: Z / sqrt(chi2_k / k)
Z = rng.standard_normal(n)
T = Z / np.sqrt(X / k)
print(f"t({k})    : var  {T.var():.4f} (theory {k/(k-2):.4f})")

And why seeding matters:

import numpy as np

a = np.random.default_rng(42).normal(size=5)
b = np.random.default_rng(42).normal(size=5)
c = np.random.default_rng(43).normal(size=5)

print("seed 42, run 1:", a.round(4))
print("seed 42, run 2:", b.round(4))
print("seed 43       :", c.round(4))
print("\nsame seed reproduces exactly:", bool(np.array_equal(a, b)))
print("different seed gives different numbers:", bool(not np.array_equal(a, c)))

Your turn

1. How would you sample from a Uniform(a, b) using a Uniform(0,1)?

2. Why can't inverse transform be applied directly to the Normal?

3. Given a stream of Uniform(0,1), how would you sample a Bernoulli(p)?

Solutions

1. Scale and shift:

X = a + (b - a)U

Via inverse transform: the CDF is F(x) = \frac{x-a}{b-a} on [a,b], so setting u = \frac{x-a}{b-a} and solving gives x = a + (b-a)u. ✓

The linear case is the one where inverse transform is completely trivial — which is why every library exposes it as a simple rescaling.

2. Because the Normal CDF

\Phi(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x e^{-t^2/2}\,dt

has no closed-form inverse. It can't even be written in elementary functions — it's the error function in disguise, and \Phi^{-1} (the probit) is only available numerically.

You can still use inverse transform with a numerical \Phi^{-1}, and scipy.stats.norm.ppf does exactly that. It's just slower and less accurate in the far tails than the alternatives.

Hence Box–Muller, which sidesteps the problem entirely by working in polar coordinates where the radial CDF is invertible; and Ziggurat, which is what production libraries use for speed.

3. Compare the uniform draw against p:

X = \begin{cases}1 & \text{if } U < p\\ 0 & \text{otherwise}\end{cases}

This works because P(U < p) = p for a Uniform(0,1) — the probability of landing in [0, p) is its length.

It's also inverse transform in the discrete case: the CDF jumps from 1-p to 1 at x = 1... or with the convention X \in \{0,1\}, from 0 to 1-p at x=0. Either way you're finding which CDF step contains U.

Everything else follows: a Binomial is a sum of these, a Geometric counts until the first success, and a categorical distribution generalises the comparison to a cumulative-sum search.

Check yourself in code

Implement inverse transform sampling for the density f(x) = 3x^2 on [0,1] and verify it against theory.

Print exactly this:

mean 0.7499
theory 0.75
P(X < 0.5) 0.125
theory 0.125

Use default_rng(1) with 500000 draws, and round every value to 4 decimal places.

import numpy as np

rng = np.random.default_rng(1)
n = 500_000

U = rng.uniform(0, 1, n)
# F(x) = x^3 on [0,1], so the inverse is u^(1/3).
X = U ** (1 / 3)

print("mean", round(X.mean(), 4))
print("theory", 0.75)

# Print the sampled P(X < 0.5) and its theoretical value F(0.5) = 0.5^3.

import numpy as np

rng = np.random.default_rng(1)
n = 500_000

U = rng.uniform(0, 1, n)
X = U ** (1 / 3)

print("mean", round(X.mean(), 4))
print("theory", 0.75)
print("P(X < 0.5)", round(float((X < 0.5).mean()), 4))
print("theory", round(0.5 ** 3, 4))

---

Everything starts from Uniform(0,1). Inverse transform sampling converts it into any distribution whose CDF you can invert, in one line; Box–Muller handles the Normal, which you can't; and many other distributions are cheapest to build directly from the relationships in §1. What remains hard is sampling from a density known only up to a constant — which is the next two lessons.

Next: how to sample from a distribution you can only evaluate, not invert.