Scientific Computing · Simulation

Rejection Sampling and the Price of the Envelope

The acceptance rate of a rejection sampler is exactly the reciprocal of its envelope constant, which is why a proposal that looks cheap can make the loop run thousands of trials per variate.

Inversion needs the quantile function, and many densities do not have one in closed form. Rejection sampling needs only the density, plus a second density that is easy to draw from and dominates it after scaling. The method is exact, requires no approximation, and returns a variate from the target distribution with probability one. The price is paid in trials: each accepted variate costs a geometrically distributed number of proposal draws and density evaluations, and the parameter of that geometric distribution is fixed by how tightly the envelope fits. A loose envelope turns an exact sampler into an expensive one, and in more than a few dimensions a loose envelope is the only kind available.

The acceptance probability is the reciprocal of the bound

Take a target density $f$ and a proposal density $q$ with $f(x) \le M q(x)$ for all $x$ and some finite $M \ge 1$. The sampler draws $Y \sim q$, draws $U$ uniform on $[0,1)$, and accepts $Y$ when $U \le f(Y) / (M q(Y))$. The acceptance probability is

$$P(\text{accept}) = \int q(y) \cdot \frac{f(y)}{M q(y)} , dy = \frac{1}{M} \int f(y) , dy = \frac{1}{M}.$$

The identity is worth reading twice. The acceptance rate does not depend on the shape of the proposal, only on the tightest constant that makes it an envelope. Improving the fit in one region while loosening it in another leaves the rate unchanged if the maximum ratio is unchanged. What matters is the single point where $f/q$ is largest.

Because trials are independent, the number of trials until acceptance is geometric with success probability $1/M$: expected trials $M$, variance $M(M-1)$, and probability of exceeding $k$ trials $(1 - 1/M)^k$. An envelope with $M = 100$ costs 100 proposal draws and 100 comparisons per accepted variate on average, and the tail is heavy enough that runs of several hundred trials occur regularly. The distribution of total work over a fixed number of accepted variates is not tightly concentrated, which is why rejection samplers are a poor fit for a fixed time budget and a good fit for a fixed sample size.

Building an envelope that actually bounds

The constant is $M = \sup_x f(x)/q(x)$, and it has to be a global supremum. A bound verified on the region where most of the mass lives is not a bound: a single point where the ratio exceeds $M$ makes the acceptance test accept more often than it should there, and the output distribution is silently wrong rather than visibly broken.

For log-concave targets, a tractable construction uses the tangent at a point. If $\log f$ is concave, its tangent at $x_0$ lies above the curve everywhere, so exponentiating the tangent gives an envelope that is proportional to a shifted exponential or, with several tangents, a piecewise exponential. Adaptive rejection sampling builds that envelope from a growing set of points, adding a point whenever a proposal is rejected, and never needs $M$ explicitly because the envelope is a normalized piecewise-exponential density whose integral is known. The construction is exact for log-concave targets and degrades to a correct but slower sampler when concavity fails, since the envelope can then fail to dominate.

A second construction, used when the target is a truncated version of a distribution with a known sampler, takes the proposal to be the untruncated density and rejects outside the truncation region. The envelope constant is then the reciprocal of the retained mass. Sampling a standard normal restricted to $[-3, 3]$ this way accepts with probability $0.9973$, which is nearly free, whereas proposing a uniform density $1/6$ on the interval gives $M = 6 \cdot 0.3989 / 0.9973 \approx 2.400$ and accepts only about $41.7%$ of trials. Both are exact; the difference is a factor of 2.4 in work per variate.

For a narrower interval the ordering reverses. Restricting the same normal to $[-1, 1]$ retains $0.6827$ of the mass, so the untruncated proposal accepts $68%$ of trials, whereas a uniform proposal on $[-1,1]$ has constant $M = 2 \cdot 0.3989 / 0.6827 \approx 1.169$ and accepts $85.5%$. The lesson is that the better proposal is whichever puts more of its mass where the target has mass, and that the comparison is arithmetic rather than a matter of taste: the ratio $f/q$ has to be evaluated at its maximum, not at the center.

The half-normal illustrates a good fit. With $f(x) = \sqrt{2/\pi}, e^{-x^2/2}$ on $x > 0$ and $q(x) = e^{-x}$, the ratio is $\sqrt{2/\pi}, e^{x - x^2/2}$, maximized at $x = 1$ where it equals $\sqrt{2e/\pi} \approx 1.316$. Acceptance is about $76%$, and the proposal costs one logarithm, the same as an inversion of the exponential.

import math, random

def half_normal():
    # Exponential envelope with rate 1; M = sqrt(2e/pi) ~ 1.316.
    M = math.sqrt(2.0 * math.e / math.pi)
    while True:
        y = -math.log(1.0 - random.random())   # Exp(1) by inversion
        if random.random() * M * math.exp(-y) <= math.sqrt(2.0 / math.pi) * math.exp(-0.5 * y * y):
            return y

The comparison in the inner loop can be rearranged to avoid the exponentials entirely: taking logarithms turns it into a quadratic test, and precomputing $\log M$ leaves one multiply and one comparison.

Squeezes and the cost of one trial

Each trial in the straightforward form evaluates the target density, which may be the expensive part. A squeeze adds a cheap lower bound $g \le f$ and accepts immediately when the uniform falls below $g(Y)/(M q(Y))$, so the target is evaluated only in the narrow band between the squeeze and the envelope. For a log-concave target with a piecewise-linear envelope, the squeeze is the chord between adjacent envelope points, and the fraction of trials that need the true density shrinks as the envelope is refined. The trade is memory and setup: more envelope points mean a tighter fit and a higher acceptance rate, at the cost of a search through the points for each trial.

The accounting that matters is total cost per accepted variate:

$$\text{cost} = M \cdot \left(c_q + c_{\text{cmp}} + p_{\text{eval}} \cdot c_f\right),$$

where $c_q$ is the proposal cost, $c_f$ the target evaluation cost, and $p_{\text{eval}}$ the probability that a trial reaches the target evaluation. A sampler with $M = 2$ and an expensive $f$ can be cheaper than one with $M = 1.2$ and a squeeze that misses often, and both can be beaten by an inversion with a table lookup if the distribution is fixed across the run.

Tails decide whether the loop terminates

The envelope condition is a statement about the tails. If the proposal has lighter tails than the target, $f/q$ is unbounded and no finite $M$ exists; the sampler will accept almost nothing in the tail and the loop will not terminate in any practical sense. A normal proposal for a Cauchy target fails this way: the ratio grows like $e^{x^2/2}/x$, so the supremum is infinite and any finite $M$ chosen by inspection is wrong at some large $x$.

In high dimensions the problem appears even when the tails are fine. Proposing a uniform on the cube $[-3,3]^d$ for a $d$-dimensional standard normal gives $M = (6/\sqrt{2\pi})^d \approx 2.394^d$, so the acceptance rate is $2.394^{-d}$: about $1/6200$ at $d = 10$ and about $1/2 \times 10^{11}$ at $d = 30$. The volume of the cube grows exponentially faster than the mass of the distribution, and rejection pays for the mismatch in trials. This is the same concentration-of-measure effect that makes naive uniform proposals useless for anything high-dimensional, and it is the reason Markov chain methods replace rejection once the dimension passes a handful.

The same exponential appears whenever the proposal is flat where the target decays, and no amount of tuning the constant removes it. A proposal that is a product of independent one-dimensional pieces has an envelope constant equal to the product of the one-dimensional constants, so a per-coordinate waste of 1.3 becomes a factor of $1.3^d$ overall. Matching the proposal to the target exactly is the only way to bring the constant back to 1, and an exact match means the proposal already is the target.

A subtler failure is a proposal that matches the target in the bulk but not in one tail, which happens with fitted envelopes. The sampler runs at a healthy acceptance rate for most draws and stalls on the rare draws from the far tail. Since the tail is exactly where the interesting rare events live, the stall is not a performance problem but a correctness one for the estimator: the samples that do arrive from the tail arrive after many rejections, and if the run is truncated by a trial budget, the tail is underrepresented.

Even a well-behaved envelope leaves the work random. For $n$ accepted variates the total trial count is negative binomial with mean $nM$ and standard deviation about $M\sqrt{n}$, so the relative fluctuation in work is $1/\sqrt{n}$ and does not shrink with a better envelope. One thread absorbs that variation by running longer. A pool of threads does not, because the slowest one sets the completion time, and the maximum of many such counts sits well above the mean. Fixed-work partitioning of a rejection sampler therefore needs over-decomposition: many more chunks than threads, so the stragglers average out.

Where rejection still wins

Rejection remains the right tool when the envelope is tight and known, when the target is available only as a density, and when exactness matters more than throughput. It is also the standard building block inside other samplers: a Metropolis-Hastings step with a symmetric proposal is a rejection step with a state-dependent envelope, and adaptive rejection sampling supplies exact draws for the conditional distributions of a Gibbs sampler when those conditionals are log-concave.

It loses to inversion whenever the quantile function is available, because inversion has a bounded cost per variate and no loop. It loses to importance sampling when the goal is an expectation rather than a sample, since importance sampling never discards a draw and instead reweights it, at the cost of a variance that depends on how well the proposal matches the integrand. The choice between the two is a choice between throwing work away and carrying a weight.

The acceptance rate also interacts with the generator. Each trial consumes a proposal variate and a uniform, so a rejection sampler advances the underlying stream by a random amount, which breaks the fixed-consumption property that makes common random numbers work across paired runs. When two configurations are compared with the same seed, a rejection step that consumes different numbers of uniforms in the two runs desynchronizes everything downstream.

References

  1. Non-Uniform Random Variate Generation — Luc Devroye; Springer
  2. Various Techniques Used in Connection with Random Digits — John von Neumann; National Bureau of Standards Applied Mathematics Series, 1951
  3. Adaptive Rejection Sampling for Gibbs Sampling — W. R. Gilks and P. Wild; Applied Statistics, 1992
  4. Monte Carlo Statistical Methods — Christian P. Robert and George Casella; Springer