Inverse-Transform Sampling and the Generalized Inverse of the CDF

A uniform generator on the half-open interval $[0,1)$ is the only primitive a random variate library really needs. Everything else is a transformation of that stream, and the oldest transformation is also the one that consumes exactly one uniform per output value: push $U$ through the inverse of the target cumulative distribution function. The method is exact in theory, cheap in the common cases, and it degrades in a specific and predictable way near the ends of the distribution. That degradation is the reason production generators carry separate tail algorithms.

The identity that makes the method work

Let $F$ be a cumulative distribution function and define the generalized inverse

$$F^{-1}(u) = \inf{x : F(x) \ge u}, \qquad 0 < u < 1.$$

For continuous $F$ this is the ordinary inverse, and the proof is one line:

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

The infimum definition matters because $F$ need not be continuous or strictly increasing. It is non-decreasing and right-continuous, so it has a left-continuous generalized inverse, and the identity survives atoms and flat regions. Choosing the infimum rather than the supremum is what makes the inequality direction come out right when $F$ has a jump exactly at $u$.

Two properties follow immediately and both are used elsewhere. The transform is monotone: if $u_1 < u_2$ then $F^{-1}(u_1) \le F^{-1}(u_2)$. And it is measure-preserving in the other direction too: $F(X)$ is uniform on $[0,1]$ whenever $X$ has continuous $F$. That second fact is the basis for goodness-of-fit tests and for the probability integral transform used to compare simulated output against a reference distribution.

Atoms become intervals

For a discrete distribution on sorted values $x_1 < x_2 < \dots$ with probabilities $p_i$, the generalized inverse is a step function. The uniform variate falls into the interval

$$\left[\sum_{j<i} p_j,; \sum_{j \le i} p_j\right)$$

with probability exactly $p_i$, so the value returned is $x_i$. This is the cumulative table method: store partial sums, draw $u$, binary search. Construction is $O(n)$ and a draw is $O(\log n)$, or $O(1)$ with the alias method, which precomputes two tables of size $n$ so that each draw needs one uniform, one multiply, one comparison and at most one extra uniform for the rejection branch.

The interval view also explains a subtlety that bites when the table is built from estimated probabilities. If the $p_i$ are rounded and do not sum to one, the last interval ends at something other than 1, and the tail of the uniform range either maps to the wrong atom or falls off the end of the table. Renormalizing after rounding fixes the sum but shifts every boundary by an amount proportional to the rounding error, so a distribution with masses near $10^{-9}$ cannot be represented faithfully in a table whose entries are stored as 32-bit floats.

import math, random

def exponential_table(lam, n=4096):
    # Equal-probability cells: boundaries are -ln(1 - i/n) / lam.
    return [-math.log1p(-i / n) / lam for i in range(n + 1)]

def exponential_draw(table, u):
    # One multiply, one floor, one table read; no logarithm per draw.
    i = int(u * (len(table) - 1))
    frac = u * (len(table) - 1) - i
    lo, hi = table[i], table[i + 1]
    return lo + frac * (hi - lo)

The interpolation inside a cell is not exact for a curved CDF; the error is second order in the cell width, which is why the cell count is chosen against the accuracy the application needs rather than against the word size.

Closed forms, and the ones that do not exist

The exponential distribution inverts in closed form. With $F(x) = 1 - e^{-\lambda x}$,

$$F^{-1}(u) = -\frac{\ln(1-u)}{\lambda}.$$

Writing $\ln(1-u)$ rather than $\ln(u)$ is deliberate. When $u$ is close to 1, the subtraction $1-u$ is exact by Sterbenz’s lemma for $u \in [0.5, 1]$, whereas $\ln(u)$ for $u$ near 1 loses relative precision because the result approaches zero. Since $1-U$ is also uniform, the reflected form costs nothing and preserves accuracy at the top of the range. The other end, $u$ near 0, is where $\ln(1-u) \approx -u$ and the formula is at its most accurate.

The normal distribution has no closed-form quantile. The practical options are a rational minimax approximation to $F^{-1}$ over the whole unit interval, accurate to roughly $10^{-15}$ in double precision, or a transformation method. Box-Muller produces two independent normals from two uniforms using one logarithm, one square root and one sine-cosine pair; the polar form (Marsaglia) replaces the trigonometric functions with a rejection loop whose acceptance probability is $\pi/4 \approx 0.785$. The inverse-transform route is slower per variate but has a property the others lack: it is monotone in $u$, so a small change in the uniform produces a small change in the output, and paired runs stay paired.

When no approximation is at hand, the quantile can be computed per draw by bracketed Newton iteration on $F(x) - u$, which needs only the CDF and the density. The bracket is what keeps it safe: Newton alone overshoots when the density is small relative to the residual, and a step outside the current bracket is replaced by bisection.

import math

def normal_cdf(x):
    return 0.5 * math.erfc(-x / math.sqrt(2.0))

def normal_quantile(u, tol=1e-15):
    lo, hi, x = -40.0, 40.0, 0.0
    for _ in range(100):
        f = normal_cdf(x) - u
        if f > 0.0:
            hi = x
        else:
            lo = x
        pdf = math.exp(-0.5 * x * x) / math.sqrt(2.0 * math.pi)
        step = f / pdf if pdf > 0.0 else 0.0
        nx = x - step
        # Newton leaves the bracket in the far tails; fall back to bisection.
        x = nx if lo < nx < hi else 0.5 * (lo + hi)
        if hi - lo <= tol * max(1.0, abs(x)):
            break
    return x

The iteration terminates on the bracket width, not on the residual, because the residual is meaningless once $F(x)$ has been rounded to the spacing of doubles near $u$. Bisection alone would need about sixty halvings to shrink the initial bracket of width 80 down to the spacing of doubles, since $80 / 2^{60}$ is below $10^{-16}$. Newton with a bisection guard reaches the same bracket in single digits for the central region and degrades toward bisection in the tails, which is the honest cost profile.

Composition and the variable number of uniforms

A mixture with weights $w_i$ and component CDFs $F_i$ has CDF $\sum_i w_i F_i$, and its inverse is the composition of two inversions: invert the discrete weights to pick a component, then invert that component. The first stage is an interval lookup, the second is whatever the component needs. Composition is how distributions with awkward shapes are handled in practice, because a finite mixture can approximate almost anything and each piece can be a distribution that already has a fast sampler.

The cost is a variable number of uniforms per output. A mixture whose components have different costs spends a different amount of the stream on each draw, and the total consumption depends on the realized weights rather than on the design. That matters when the stream is being shared between two runs that must stay synchronized, and it matters for reproducibility: if the component choice consumes a uniform and the component draw consumes another, a change to the component ordering shifts every subsequent value in the stream.

The tail behaves the way the heaviest component dictates. A mixture that puts weight $10^{-3}$ on a heavy-tailed component still has that component’s tail, and inversion of the mixture inherits the truncation problem described below at exactly the point where the heavy component’s mass becomes smaller than the spacing of doubles. Increasing the table resolution does not help, because the limit is the resolution of the uniform, not the resolution of the table.

Where the tails run out of bits

A double-precision uniform has 53 bits of mantissa, so the spacing of representable values near 1 is $2^{-53} \approx 1.1 \times 10^{-16}$. The largest value below 1 is $1 - 2^{-53}$. Any quantile beyond that point is unreachable by inversion: the event $U > 1 - 2^{-53}$ has probability about $1.1 \times 10^{-16}$, and asking for the $1 - 10^{-18}$ quantile of a normal returns the same value as the $1 - 10^{-16}$ quantile. The distribution is truncated, silently, at roughly 8.2 standard deviations.

For a simulation that estimates a mean, truncation at eight sigma is irrelevant; the omitted mass contributes on the order of $10^{-16}$ to any expectation. For a rare-event study, a tail-risk calculation, or a likelihood that is dominated by its extreme observations, it is fatal. The standard remedy is to sample the tail directly: draw the exceedance probability from a uniform on the small interval that remains, then invert the tail expansion of the distribution, which for the normal is the asymptotic series for the Mills ratio. The same reasoning applies to the discrete case, where a cumulative table in double precision cannot represent masses below $2^{-53}$ at all.

A related trap is the open versus closed interval. Some generators return values in $[0,1)$, some in $(0,1)$, and some can return exactly 0 or 1 depending on how the bits are assembled. A generator that can return 0 makes $\ln(u)$ infinite; a generator that can return 1 makes the normal quantile infinite. Guarding with a resample costs a branch and destroys the fixed one-uniform-per-variate accounting that makes inversion attractive for antithetic variates and common random numbers, where each pair must consume the stream in lockstep.

Why inversion is still the default

Rejection methods can be faster per variate for some distributions, and composition methods beat inversion for mixtures. Inversion wins on three counts that matter in simulation rather than in a single draw. It consumes exactly one uniform, so streams stay aligned across runs. It is monotone, so negatively correlated inputs produce negatively correlated outputs, which is the mechanism behind variance reduction by antithetics. And it is the only method that makes the variance of the estimator depend on the variate through $F$ alone, which allows stratification: splitting the unit interval into $k$ equal strata and drawing once from each gives a sample whose empirical distribution is closer to $F$ than $k$ independent draws, without changing the per-draw cost.

The cost is that every distribution needs its own inverse, and the ones without a closed form need either a fitted approximation with documented error bounds or a numerical solve per draw. The numerical solve is the fallback that always works: bracket the quantile by bisection on $F$, then refine with Newton using the density, and stop when the bracket is at the resolution of the output type. It costs a handful of CDF evaluations per variate, which is the price of not having an approximation whose error has to be trusted.

The next question is where the uniforms themselves come from, since the whole construction inherits their defects. A generator with lattice structure produces correlated tuples, and inversion maps that correlation into the output without removing it; the lattice structure of linear congruential generators is the canonical case.