How to draw a sample

You have uniforms; you need this distribution

Imagine you are handed a black box that returned 0.3719, then 0.0288, then 0.9144, each an independent uniform draw on the unit interval. That is all the randomness a computer can honestly manufacture. Now suppose the task is to simulate the waiting time between two radioactive decays, the height of a person drawn from a population, or the next token in a sentence. Those are not uniform, and the whole usefulness of simulation depends on converting one into the other without bias. The sample you produce must have the right distribution, not merely the right range.

The classical framing is this. You are given a target distribution F — maybe only through its density f up to an unknown constant, which is the normal state of affairs in Bayesian inference. You are given an unlimited supply of independent uniform variables, and often a second distribution g that you can sample. The question is how to turn draws from the uniform, or from g, into draws from F. Different answers give different costs. Inversion is exact and uses one uniform per sample but needs the quantile function. Rejection is exact and needs almost nothing but wastes draws. Importance sampling does not even return a sample from F; it returns a weighted collection from which expectations under F can still be estimated.

It is worth saying why we care about exactness at all. A simulator that draws from slightly the wrong distribution produces slightly wrong answers, and the error is invisible because it never averages out — more samples shrink the Monte Carlo noise around the wrong value. In a filter tracking a robot, a biased sampler means a drifting estimate; in a reinforcement-learning gradient, it means the policy improves toward the wrong objective. The three methods below are the standard ways to keep the distribution honest while keeping the arithmetic simple.

The plan for this part follows the three bridges. First inverse-CDF sampling, the cleanest of the three, with a live conversion of uniforms into samples from a normal, an exponential or a beta, checked against the target density. Then rejection sampling with an envelope you can drag, an acceptance meter, and a visible pile of rejected points. Then importance sampling, where the object of study is the weight histogram and the effective sample size, and where a deliberately bad proposal makes the method collapse on command. The last sections collect the formulas and send you on to Monte Carlo integration, which consumes exactly these samplers.

Push a uniform through the quantile function

Start with the target CDF $F(x)=P(X\le x)$. It is non-decreasing and runs from zero to one, so it can be inverted: for a number u in the unit interval, define F^{-1}(u) to be the smallest x with $F(x)\ge u$. Feed a uniform variable U into that inverse and you get a variable with CDF F. The proof is one line. Because F is monotone, the event $F^{-1}(U)\le x$ is the same as the event $U\le F(x)$; and since U is uniform on [0,1], the probability of $U\le F(x)$ is exactly F(x). That is what it means for the new variable to have distribution F.

This is the probability integral transform, and it is a two-way street. Run the argument backwards and a variable with continuous CDF F, pushed through F itself, becomes uniform — the fact that makes probability plots and copulas work. The forward direction is what we use for sampling. The picture in the demo is the geometry: a uniform on the vertical axis, a horizontal line to the CDF curve, and a vertical drop to the horizontal axis. Uniformly spaced levels on the y-axis, once mapped through the curve, land densely where the CDF rises steeply, which is exactly where the density is large.

What does it cost? One uniform per sample, no rejections, no tuning. The catch is the quantile function. For a normal distribution there is no elementary closed form for F^{-1}, so a library either inverts the CDF numerically — a bisection or Newton search, as the Prob distributions here do — or falls back on a special-purpose construction such as Box–Muller. For the exponential, $F^{-1}(u)=-\log(1-u)/\lambda$, and for the uniform it is a linear map, both essentially free. Discrete targets are even easier: ordering the outcomes and comparing the uniform against the running cumulative probabilities is inversion, and it is precisely the algorithm a language model uses to sample a token from a softmax.

Two practical cautions. First, the inversion must be numerically stable near the ends: F^{-1} blows up or saturates as u approaches zero or one, and a quantile routine that bisects a CDF with flat tails can return garbage if the bracket is wrong. Second, inversion delivers independent samples only because the uniforms were independent. Push a single uniform through many different quantile functions and you get a perfectly dependent vector — the basis of copulas, but a trap if you expected a random sample. In the demo below, switch between the three families and draw: the histogram of the mapped uniforms sits under the target density, because that is what the theorem promises.

Notice that the two panels tell the same story from opposite ends. The upper panel is the algorithm, mechanical and deterministic once the uniforms are chosen; the lower panel is the certificate, statistical and noisy. If the quantile function were wrong by even a little, the upper panel would look fine and the lower would be tilted, too heavy in one tail or shifted off centre. That is why simulation code is tested against moments and quantiles rather than eyeballed: a sampler is correct exactly when its output density matches the target.

Enclose the target, throw darts, keep what lands inside

Inversion needs the quantile function. Rejection sampling needs only the density, and it needs it only up to a constant — no integral, no inversion, no normalising constant. Pick a proposal distribution g that you already know how to sample and that is positive wherever the target f is. Choose a constant M large enough that $M\,g(x)\ge f(x)$ for every x. That curve $M\,g$ is the envelope: the proposal, scaled up until it sits above the target everywhere. Then repeat the same two steps. Draw $X\sim g$ and an independent $U\sim\mathrm{Uniform}(0,1)$; accept X if $U\le f(X)/\bigl(M\,g(X)\bigr)$, otherwise throw it away.

Why is the accepted point distributed as f? A proposed point at location x is accepted with probability f(x)/(M g(x)), and the chance of proposing near x is proportional to g(x). Multiply the two: the density of accepted points is proportional to $g(x)\cdot f(x)/(M g(x))=f(x)/M$. The proposal cancels, leaving the target. Normalising over all x gives f exactly. The same computation gives the acceptance probability: integrating the acceptance probability against g yields $\frac{1}{M}\int f=\frac{1}{M}$, so the method keeps an average of one draw in M.

The envelope constant is therefore the entire efficiency budget. Make M as small as the constraint allows — the tightest envelope is $M^*=\sup_x f(x)/g(x)$ — and the rejection rate is 1-1/M^*. A proposal shaped like the target wastes almost nothing; a proposal that is far too wide, or centred in the wrong place, wastes almost everything. Rejection sampling is also unforgiving about the support: if g is zero somewhere f is not, that region of the target is unreachable and the result is silently biased. The demo flags this directly — an envelope that dips below the target draws a warning, because then the accept rule no longer produces the target.

There is a hard limitation hiding in the constant. In one dimension, fitting an envelope is easy. In ten dimensions, the volume between a decent proposal and a target is enormous, and M^* grows roughly exponentially with the dimension; the acceptance probability collapses to nothing. This is the curse of dimensionality, and it is the reason rejection sampling is used in low dimensions and replaced by Markov chain Monte Carlo — the subject of a later part — in high ones. A related cousin, adaptive rejection sampling, rebuilds the envelope on the fly when f is log-concave, squeezing M down automatically.

Drag the envelope handle on the canvas to scale M, and use the sliders to shift and widen the Gaussian proposal. Watch three things move together: the envelope height, the fraction of points that land under the target (accepted, in blue) rather than in the gap above it (rejected, in grey), and the resulting histogram. Tighten the envelope until it just grazes the target and the acceptance meter climbs toward its theoretical ceiling of 1/M. Pull it up into a tall, flat lid and almost every draw is wasted.

Rejection sampling is the honest workhorse of low-dimensional simulation, and its logic reappears elsewhere under other names. Deciding whether a proposal is acceptable is a coin flip whose bias is the likelihood ratio, and the same ratio governs acceptance in Metropolis–Hastings. The practical skill it teaches is the one that carries over: to sample something hard, sample something easy and correct the answer for the difference. When the correction is a hard accept/reject decision, you get this method; when it is a weight, you get the next one.

Keep every draw, but weight it

Rejection throws away the points it does not want. Importance sampling keeps them all and records how much each one should count. The change is small in code and large in effect: instead of producing a sample from f, it produces a sample from a proposal g together with a weight w(x)=f(x)/g(x), and every expectation is then taken as a weighted average. The identity that makes it work is a one-line change of measure. For any function h whose expectation under f exists, write the integral against f, multiply and divide by g, and read the result as an integral against g.

So draw $x_1,\dots,x_n\sim g$ and estimate $\mathbb{E}_f[h]$ by the average of h(x_i)w(x_i). That estimator is unbiased, provided g is positive wherever h f is not zero — the same support condition as rejection, with the same silent bias if violated. Its variance is $\frac{1}{n}\bigl(\mathbb{E}_g[h^2w^2]-\mathbb{E}_f[h]^2\bigr)$, and now the proposal has real leverage. The variance is small when g is large exactly where |h|f is large. In the ideal case $g^*\propto|h|f$ the variance is zero: every weight is identical and the estimate is exact before any noise.

In practice the target density is usually known only up to a constant, so the plain average cannot be formed and one uses the self-normalised estimator $\widehat I_{\mathrm{SN}}=\sum_i w_i h(x_i)/\sum_i w_i$. It is slightly biased for finite n but consistent, never needs the normalising constants, and is what the demo computes. The quality of that estimator is governed not by n but by how evenly the weight is spread. Define the effective sample size as $\mathrm{ESS}=(\sum_i w_i)^2/\sum_i w_i^2$, the number of equally weighted draws that would carry the same information. When the weights are all equal, $\mathrm{ESS}=n$; when one draw holds all the weight, $\mathrm{ESS}=1$, no matter how many samples you drew.

Degeneracy is the characteristic failure. If g rarely lands where f has mass, then a few lucky draws fall in the target's core, acquire enormous weights, and dominate the average; the rest contribute nothing. The estimate then depends on whether those lucky draws happened, which is a variance problem that no amount of extra sampling fixes — adding more draws from the same bad proposal adds more negligible-weight points, not more information. The diagnostic is exactly the weight histogram below: a healthy run puts almost all weight near the same value, while a bad proposal produces a long tail with a spike near zero and a handful of giants.

The demo samples from a Gaussian proposal and estimates the mean of a standard normal target using h(x)=x. Move the proposal centre away from the target or shrink its width and the weight cloud stretches out, the largest weight swallows the total, and the ESS meter drains even though the sample count never changes. Move it back and the weights flatten. This is the picture to keep in mind whenever you see an importance weight: it is a correction factor, and corrections that vary by orders of magnitude are a warning that the two distributions barely overlap.

Importance sampling is not a sampler in the sense of the first two sections — it never yields an unweighted draw from f. What it yields is a consistent estimator of any expectation under f, and that is often all that was wanted. The weighted cloud of particles is the seed of sequential importance sampling and the particle filter, where the proposal is built one observation at a time and the weights are multiplied by each likelihood. The ratio f/g is also the likelihood ratio that appears in off-policy policy gradients and in the clipped objective of proximal policy optimisation, where it is the same correction repurposed as a training signal.

Every simulator is one of these three in disguise

Every formula in one place

Where to go deeper

• Sheldon Ross, Simulation, 5th edition — the clearest textbook treatment of the inverse-transform and acceptance-rejection methods, with the algorithms written out fully.
• Jun S. Liu, Monte Carlo Strategies in Scientific Computing, 2001 — the standard reference on importance sampling, sequential importance sampling and the theory of effective sample size.
• Art B. Owen, Monte Carlo Theory, Methods and Examples, 2013 — an accessible and rigorous chapter on importance sampling with variance analysis and the ESS diagnostic.
• Christian P. Robert and George Casella, Monte Carlo Statistical Methods, 2nd edition — rejection and importance sampling placed in the larger arc toward MCMC.
• Reuven Y. Rubinstein and Dirk P. Kroese, Simulation and the Monte Carlo Method, 3rd edition — the crossover and change-of-measure view of importance sampling, useful when the proposal is chosen for variance reduction.
• The volume glossary — a routing table from each of these ideas into the AI, vision and robotics pages that depend on them.