65. Rejection sampling

Inverse transform needs F^{-1}. Often you have only f — or worse, something proportional to f with an unknown normalising constant.

Rejection sampling needs neither. It only needs to evaluate the density, and it works for any distribution you can bound.

The algorithm

Pick a proposal distribution q you can sample from, and a constant M with

f(x) \le M\,q(x) \quad \text{for all } x

Then repeat:

1. Draw X \sim q.
2. Draw U \sim \text{Uniform}(0,1).
3. Accept X if U \le \dfrac{f(X)}{M\,q(X)}; otherwise discard and repeat.

Accepted values are exact draws from f. Not approximate — exact.

Why it works

Picture the 2-D region under M q(x). Step 1 picks a horizontal position from q; step 2 picks a height uniformly within the envelope at that position. Together, the pair (X, U \cdot Mq(X)) is uniform over the area under Mq.

Accepting exactly when the point falls under f keeps a uniform sample from the area under f — and the x-coordinates of points uniform under a density are distributed according to that density.

Formally:

P(X \le x \mid \text{accept}) = \frac{\int_{-\infty}^x q(t)\frac{f(t)}{Mq(t)}dt}{\int_{-\infty}^\infty q(t)\frac{f(t)}{Mq(t)}dt} = \frac{\frac{1}{M}\int_{-\infty}^x f}{\frac{1}{M}\int f} = F(x)

The q cancels completely, which is the whole trick.

The acceptance rate

P(\text{accept}) = \frac{1}{M}

and the expected number of proposals per accepted sample is M — a Geometric waiting time (§1).

So you want M as small as possible, which means q should hug f as tightly as you can manage. M must be at least 1, with equality only if q = f.

The curse of dimensionality

Rejection sampling is excellent in one or two dimensions and collapses in high ones.

The reason: M typically grows exponentially with dimension. If the proposal mismatches the target by a factor c per dimension, then in d dimensions M \sim c^d and the acceptance rate is c^{-d}.

A method accepting 50% in one dimension accepts 0.5^{20} \approx 10^{-6} in twenty — a million proposals per sample.

This is why MCMC exists. Markov chain methods (§8) explore the target locally rather than proposing globally, and their cost scales far better with dimension. Rejection sampling remains valuable as a building block — inside Gibbs samplers, or for low-dimensional sub-problems.

Unnormalised targets

Here is the property that makes it genuinely useful.

If you only know \tilde f \propto f, the method still works: find M with \tilde f(x) \le Mq(x) and accept with probability \frac{\tilde f(X)}{Mq(X)}.

The normalising constant never appears. It cancels in the ratio, exactly as q did.

That's precisely the Bayesian situation from §7: the posterior is known up to the intractable evidence term p(\text{data}), and rejection sampling doesn't need it.

Worked example

Sample from f(x) = 2x on [0,1] using a uniform proposal.

Proposal: q(x) = 1 on [0,1].

Bound: \max f = f(1) = 2, so M = 2 works.

Algorithm: draw X, U uniform; accept X if U \le \frac{2X}{2 \cdot 1} = X.

The acceptance condition simplifies beautifully to U \le X.

Acceptance rate: 1/M = 0.5 — half the proposals are kept, so about 2 proposals per sample.

Check: E[X] = \int_0^1 x(2x)dx = 2/3 \approx 0.667, and the sampler gives 0.667.

Now with a bad proposal. Use $q = $ Beta(0.5, 0.5), which piles mass at both endpoints while f is largest at 1 and zero at 0. The required M jumps, and the acceptance rate falls sharply — same target, same correctness, far more work. The proposal doesn't affect correctness, only cost.

Doing it in Python

The f(x) = 2x example:

import numpy as np

rng = np.random.default_rng(0)

def rejection_sample(n, rng):
    """Sample f(x) = 2x on [0,1] with a uniform proposal, M = 2."""
    accepted, proposals = [], 0
    while len(accepted) < n:
        batch = 2 * (n - len(accepted)) + 100
        X = rng.uniform(0, 1, batch)
        U = rng.uniform(0, 1, batch)
        keep = X[U <= X]                     # U <= f(X)/(M q(X)) = X
        proposals += batch
        accepted.extend(keep.tolist())
    return np.array(accepted[:n]), proposals

X, proposals = rejection_sample(200_000, rng)

print(f"mean            : {X.mean():.4f}   theory {2/3:.4f}")
print(f"P(X < 0.5)      : {(X < 0.5).mean():.4f}   theory {0.25:.4f}")
print(f"acceptance rate : {200_000/proposals:.4f}   theory {1/2:.4f}")

An unnormalised target — the Bayesian case:

import numpy as np
from scipy.stats import beta

rng = np.random.default_rng(1)

# Posterior for a coin: Beta(3,2) prior x 6 heads in 9 flips, UNNORMALISED
def unnormalised(p):
    return p**2 * (1 - p)**1 * p**6 * (1 - p)**3      # prior x likelihood

grid = np.linspace(0, 1, 10_001)
M = unnormalised(grid).max() * 1.001                   # bound with a uniform proposal

n, accepted, proposals = 200_000, [], 0
while len(accepted) < n:
    batch = 5 * (n - len(accepted)) + 100
    X = rng.uniform(0, 1, batch)
    U = rng.uniform(0, 1, batch)
    accepted.extend(X[U <= unnormalised(X) / M].tolist())
    proposals += batch
X = np.array(accepted[:n])

exact = beta(3 + 6, 2 + 3)      # the conjugate answer, for checking
print(f"sampled mean : {X.mean():.4f}   exact Beta(9,5) mean {exact.mean():.4f}")
print(f"sampled var  : {X.var():.4f}   exact var  {exact.var():.4f}")
print(f"acceptance   : {n/proposals:.4f}")
print("\nThe normalising constant was never computed -- it cancels in the ratio.")

How the proposal choice affects cost, never correctness:

import numpy as np
from scipy.stats import norm, cauchy, uniform

rng = np.random.default_rng(2)

# Target: standard normal, restricted to [-5, 5] so a uniform proposal works
target = lambda x: norm.pdf(x)

proposals = {
    "uniform on [-5,5]": (lambda k: rng.uniform(-5, 5, k), lambda x: uniform.pdf(x, -5, 10)),
    "Cauchy(0,1)":       (lambda k: rng.standard_cauchy(k), lambda x: cauchy.pdf(x)),
    "Normal(0, 1.5)":    (lambda k: rng.normal(0, 1.5, k), lambda x: norm.pdf(x, 0, 1.5)),
}

grid = np.linspace(-5, 5, 20_001)
for name, (sampler, density) in proposals.items():
    M = (target(grid) / density(grid)).max() * 1.001
    X = sampler(300_000)
    X = X[(X > -5) & (X < 5)]
    U = rng.uniform(0, 1, len(X))
    kept = X[U <= target(X) / (M * density(X))]
    print(f"{name:>20}: M = {M:>6.3f}, acceptance {len(kept)/len(X):>7.4f}, "
          f"sample mean {kept.mean():+.4f}")

print("\nEvery proposal gives the right answer. They differ only in how many")
print("draws get thrown away.")

The dimensionality collapse:

import numpy as np

rng = np.random.default_rng(3)

# Target: uniform on the unit ball. Proposal: uniform on the cube [-1,1]^d.
# Acceptance = volume ratio of ball to cube.
print(f"{'dimension':>10} {'acceptance rate':>18} {'proposals per sample':>22}")
for d in (1, 2, 3, 5, 10, 20):
    n = 200_000
    X = rng.uniform(-1, 1, size=(n, d))
    inside = (np.sum(X**2, axis=1) <= 1).mean()
    per = f"{1/inside:,.0f}" if inside > 0 else "> 200,000"
    print(f"{d:>10} {inside:>18.6f} {per:>22}")

print("\nBy d = 20 the ball occupies essentially none of the cube.")
print("This is why MCMC replaces rejection sampling in high dimensions.")

Your turn

1. For f(x) = 2x on [0,1] with uniform proposal, what's the acceptance rate?

2. Why does rejection sampling work with an unnormalised target?

3. Your acceptance rate is 0.1%. What should you change?

Solutions

1. 50%.

The acceptance rate is 1/M, and we need f(x) \le Mq(x) with q = 1 on [0,1]. Since \max_{[0,1]} 2x = 2, the smallest valid M is 2, giving

P(\text{accept}) = \frac{1}{2} = 0.5

Geometrically: the envelope is the rectangle [0,1] \times [0,2] with area 2, and the region under f is a triangle of area 1. The ratio is 1/2.

On average you need 2 proposals per accepted sample.

2. Because the constant cancels in the acceptance ratio.

Write f(x) = \tilde f(x)/Z where Z = \int \tilde f is unknown. Choose M with \tilde f(x) \le Mq(x) and accept with probability \frac{\tilde f(X)}{Mq(X)}.

The proof of correctness normalises the accepted distribution by dividing by the total acceptance probability, and Z appears identically in numerator and denominator:

P(X \le x \mid \text{accept}) = \frac{\int_{-\infty}^x \tilde f}{\int_{-\infty}^{\infty}\tilde f} = F(x)

This is the property that makes it useful for Bayesian inference. The posterior p(\theta \mid \text{data}) \propto p(\text{data} \mid \theta)p(\theta) has an intractable normalising constant (§7), and rejection sampling never needs it.

The one cost: without Z you can't easily know how good your M is in absolute terms, so the acceptance rate must be measured empirically.

3. A 0.1% rate means 1,000 proposals per sample — the proposal is a poor match for the target.

Options, roughly in order:

• Choose a better q. The most direct fix. Match the target's location, scale, and especially its tails — a proposal with lighter tails than the target forces a huge M (or makes a valid M impossible).
• Check M isn't needlessly large. You want the smallest valid bound, M = \sup_x f(x)/q(x). A lazily chosen M wastes proposals for no reason.
• Use adaptive rejection sampling if the log-density is concave — it builds a piecewise envelope that tightens as it goes.
• Switch method. If the problem is high-dimensional, no proposal will save you: M grows exponentially with d. Move to MCMC, which explores locally, or importance sampling (next lesson), which weights rather than discards.

That last point is the practical conclusion of this lesson. A 0.1% acceptance rate in a 10-dimensional problem is not a tuning failure — it's the method telling you it's the wrong tool.

Check yourself in code

Implement rejection sampling for f(x) = 2x on [0,1] and verify both the distribution and the acceptance rate.

Print exactly this:

mean 0.6667
theory 0.6667
acceptance 0.4993
theory 0.5

Use default_rng(0), generate 200000 accepted samples, propose from Uniform(0,1) with M = 2, and round every value to 4 decimal places.

import numpy as np

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

accepted, proposals = [], 0
while len(accepted) < n:
    batch = 2 * (n - len(accepted)) + 100
    X = rng.uniform(0, 1, batch)
    U = rng.uniform(0, 1, batch)
    # Accept when U <= f(X) / (M q(X)) = 2X / 2 = X
    accepted.extend(X[U <= X].tolist())
    proposals += batch

X = np.array(accepted[:n])
print("mean", round(X.mean(), 4))

# Print the theoretical mean 2/3, the empirical acceptance rate n/proposals,
# and its theoretical value 1/M.

import numpy as np

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

accepted, proposals = [], 0
while len(accepted) < n:
    batch = 2 * (n - len(accepted)) + 100
    X = rng.uniform(0, 1, batch)
    U = rng.uniform(0, 1, batch)
    accepted.extend(X[U <= X].tolist())
    proposals += batch

X = np.array(accepted[:n])
print("mean", round(X.mean(), 4))
print("theory", round(2 / 3, 4))
print("acceptance", round(n / proposals, 4))
print("theory", 0.5)

---

Rejection sampling draws from an envelope and keeps the points that fall under the target, giving exact samples from any density you can evaluate and bound — including unnormalised ones, which is why it suits Bayesian posteriors. Its cost is M proposals per sample, and since M grows exponentially with dimension, it belongs to low-dimensional problems.

Next: using samples to compute integrals.