Inverse & Rejection Sampling theory

The Inverse Transform method (Probability Integral Transform) states that if X is a continuous random variable with strictly increasing cumulative distribution function F(x), then U = F(X) is uniformly distributed on [0, 1]. Inversely, if U ~ Uniform(0, 1), then X = F^(-1)(U) has exact cumulative distribution function F. The proof is immediate: P(X ≤ x) = P(F^(-1)(U) ≤ x) = P(U ≤ F(x)) = F(x). The method is universally exact and achieves 100% computational efficiency with zero rejected draws, but requires that F^(-1) exist in closed analytical or rapidly computable numerical form.

Introduced by John von Neumann in 1951, Acceptance-Rejection sampling generates variates from a target distribution with density f(x) using an accessible proposal distribution with density g(x). A bounding constant c is chosen such that f(x) ≤ c g(x) for all x in the support, so that c g(x) forms a majorizing envelope over f(x). The algorithm proceeds in two steps: (1) Generate candidate Y ~ g(x) and independent threshold U ~ Uniform(0, 1); (2) If U ≤ f(Y) / [c g(Y)], accept X = Y; otherwise, reject Y and repeat. The accepted variates follow density f(x) exactly.

Let A denote the event that a candidate proposal is accepted: A = {U ≤ f(Y) / (c g(Y))}. Conditioning on the candidate value Y = y yields P(A | Y = y) = f(y) / (c g(y)). Integrating over the proposal distribution gives P(A) = ∫ P(A | Y = y) g(y) dy = ∫ [f(y) / (c g(y))] g(y) dy = (1/c) ∫ f(y) dy = 1/c, since f(y) is a valid probability density integrating to 1. The number of proposals required to obtain one accepted variate follows a Geometric distribution with parameter p = 1/c, having expected value E[N] = 1/p = c.

To maximize sampling efficiency (acceptance rate 1/c), the bounding constant c must be minimized subject to the envelope condition: c_opt = sup_{x} [f(x) / g(x)]. If c < c_opt, the envelope fails to majorize f(x), distorting the target distribution and producing biased samples. If c is chosen unnecessarily large (e.g., in high-dimensional spaces), the acceptance probability 1/c drops exponentially toward zero, rendering the algorithm computationally infeasible (the curse of dimensionality).

Developed by George Marsaglia and Wai Wan Tsang in 2000, the Ziggurat algorithm covers the target probability density (such as the standard Gaussian curve) with a stack of horizontal rectangular layers (resembling a ziggurat). The layers are constructed to have identical areas. For 98% of proposals falling within the completely enclosed core of a rectangle, the variate is accepted immediately without evaluating expensive exponential or transcendental functions, making it the fastest known algorithm for normal and exponential generation.

Introduced by Gilks and Wild in 1992, Adaptive Rejection Sampling constructs piecewise linear upper and lower bounding envelopes in the log-density domain ln f(x). Because f(x) is log-concave (d^2 ln f(x) / dx^2 ≤ 0), tangents to the log-density curve lie strictly above the curve, forming an envelope whose exponent is piecewise exponential. Whenever a candidate is rejected, the evaluated point is added to update the envelope, causing the envelope to converge progressively to f(x) and driving acceptance efficiency toward 100%.

Classical rejection sampling requires constructing a global majorizing envelope c g(x). In d dimensions, the volume of the target density relative to any standard proposal envelope typically shrinks as (1/c)^d, causing acceptance rates to collapse to near zero. MCMC algorithms (such as Metropolis-Hastings) bypass the need for a global envelope by generating a correlated random walk that explores the target distribution locally based only on local density ratios f(x')/f(x), making simulation viable in high-dimensional Bayesian models.