A language model ends with 50,000 logits and has to pick one. Sampling is the arithmetic of that choice — and the same recipes estimate integrals, train VAEs, and walk Bayesian posteriors. Every one of them starts from U ~ Uniform(0, 1).
The CDF turns a value into a probability. Running it backwards — pick a random height on the y-axis, read across to the x-axis — turns a probability into a value. X = F⁻¹(U) is an exact sample, and the steep parts of the curve catch most of the random heights.
X = F⁻¹(U)Pick a simple proposal q that is easy to sample, and inflate it by M until the curve M·q sits above the target p everywhere. Draw x ~ q; accept it with probability p(x)/(M·q(x)). Accepted points are exact samples from p, but 1 − 1/M of the work is thrown away.
accept if u < p(x)/(M·q(x))When you need E_p[f] but can only draw from q, multiply each sample by w = p(x)/q(x). Samples the target values more count more; samples it dislikes count less. Or skip reweighting entirely and walk a Markov chain whose long-run visits are the target — that is Metropolis–Hastings.
E_p[f] ≈ (1/N)·Σ f(xᵢ)·w(xᵢ)A computer makes exactly one cheap flavor of randomness — U ~ Uniform(0, 1) — so every sampler is a recipe for turning that uniform draw into the distribution you actually want: transform it (inverse CDF), filter it (rejection), reweight it (importance sampling), or walk a chain until it forgets where it started (MCMC).
By the end you will be able to sample any distribution from uniforms, predict acceptance rates, correct an expectation with importance weights, run and diagnose a Metropolis–Hastings chain, and tune temperature, top-k and top-p on a language model — plus see why the reparameterization trick is what makes VAEs trainable.
A language model finishes with 50,000 logits and has to pick one. That pick is sampling — and the same handful of recipes shows up in training, estimation and exploration.
Sampling has four jobs in modern AI. Generation: language models, diffusion models and GANs produce output by drawing from a distribution. Training: stochastic gradient descent samples mini-batches, dropout samples which neurons to switch off, and data augmentation samples random transformations. Estimation: quantities with no closed form — the expected loss over a data distribution, the evidence in Bayesian inference — get approximated by averaging samples. Exploration: MCMC walks around a posterior, and Thompson sampling picks between options with unknown payoffs.
All of it starts from one primitive: a uniform random number U ~ Uniform(0, 1). Every sub-interval of equal length is equally likely, its mean is ½ and its variance is 1/12. To pick one of n items, take ⌊n·U⌋; to land anywhere in [a, b] with equal density, take a + (b − a)·U. The rest of this lesson is recipes that turn that single uniform draw into draws from whatever distribution you actually need.
Each experiment draws n numbers from a standard normal and estimates the variance two ways. The biased estimator lands on the wrong value forever; both get tighter as n grows.
This is why (n − 1) appears in the sample variance: dividing by n systematically undershoots, because the sample mean is fitted to the same data. The undershoot is exactly a factor (n − 1)/n.
Plain English: an average of many independent draws has the same mean as one draw, but a much smaller spread — shrinking with the square root of how many draws you average.
Notice what is not in that formula: the dimension of x. A grid with 10 points per axis covers 10² cells in 2-D but 10¹⁰⁰ cells in 100-D. Monte Carlo error ignores all of that, which is why it takes over exactly where grids die.
There is a second lesson in the lab above: an estimator can be wrong in two different ways. Bias is a mean that misses the target; variance is a spread around its own mean. They combine into the mean squared error:
You need a Monte Carlo estimate to be 10× more precise. How many more samples do you need?
If you can write down the cumulative distribution function and invert it, one uniform draw becomes an exact sample from any distribution — no rejection, no approximation, no wasted randomness.
The cumulative distribution function (CDF) is the running total of probability: F(x) = P(X ≤ x). It starts at 0, climbs without ever going down, and ends at 1. The inverse-CDF method — inverse transform sampling — says: draw u ~ Uniform(0, 1) and return x = F⁻¹(u). The output is distributed exactly as F.
Plain English: pick a uniformly random height on the y-axis of the CDF, then walk horizontally to the curve and drop down to the x-axis. Where the CDF is steep, many heights land in the same short stretch of x; where it is flat, very few do. The steepness of the CDF is the density — so popular values collect the most samples.
One random height u on the CDF’s y-axis, read across to the x-axis. Steep parts of the CDF catch most of the heights, so they collect most of the samples.
The formula for the exponential is x = −ln(1 − u)/λ. The triangular case is even simpler: x = √u. The Cauchy is x = tan(π(u − ½)) — its heavy tails are real, and its running mean never settles because the expected value does not exist.
The proof is one line of bookkeeping: applying the non-decreasing F to both sides of an inequality preserves it, and a uniform variable lands below a number with probability equal to that number. The interesting part is the arithmetic for each distribution.
F(x) = 1 − e^(−λx). Set F(x) = u and solve: e^(−λx) = 1 − u, so x = −ln(1 − u)/λ. Since 1 − U has the same uniform distribution as U, x = −ln(u)/λ is equally valid.2x on [0, 1] has F(x) = x², so the inverse is x = √u. Small uniforms map to even smaller values — exactly the rising density.u = 0.62 lands in the second category, because 0.5 < 0.62 ≤ 0.8. This is exactly how sample_categorical works.F(x) = ½ + arctan(x)/π inverts to x = tan(π(u − ½)). The tails are so heavy that the mean does not exist: u = 0.975 already gives x ≈ 12.7.The normal distribution has no closed-form inverse CDF, which is why Box–Muller builds it from two uniforms instead: z = √(−2·ln u₁)·cos(2πu₂) is a standard normal sample. With u₁ = 0.5, u₂ = 0 it gives z = 1.177. Two uniforms, two independent normals — no rejection, no tables.
An exponential has λ = 2. With u = 0.5, what is x = −ln(1 − u)/λ?
When you can evaluate the target density but cannot invert its CDF, build a ceiling above it with a distribution you can sample, then accept or reject each proposal.
Rejection sampling needs three ingredients: a target p(x) you can evaluate (even up to a constant), a proposal q(x) you can sample from, and an envelope constant M big enough that M·q(x) sits above p(x) everywhere. Draw x ~ q(x), draw u ~ Uniform(0, 1), and accept x when u < p(x)/(M·q(x)). Otherwise throw the point away and try again.
Plain English: sprinkle points uniformly under the ceiling, keep the ones that land under the target curve, and pretend the rest never existed. The kept points are exact samples from p — rejection sampling is not approximate. It is just wasteful, and it gets exponentially more wasteful in higher dimensions.
Proposals fall from a uniform distribution. A point is kept only if it lands under the target curve — green accepted, red rejected. Tune the envelope and watch the acceptance rate.
Acceptance is exactly the area under p divided by the area under M·q — which is 1/M. It falls exponentially with dimension, so beyond a few dimensions this exact method is replaced by MCMC.
That last line is why rejection sampling is a 1-D-to-3-D tool. High-dimensional posteriors need a method that does not require a ceiling over the whole space — which is what the MCMC chapter builds.
One more numeric habit: M must clear the peak everywhere, not just on average. For Beta(2, 5), choosing M = 2.0 leaves the top of the curve poking through the ceiling. Wherever p(x) > M·q(x), the acceptance probability exceeds 1 and those proposals are always kept — the sampler silently over-samples the peak and the mean drifts away from 2/7 ≈ 0.286. The lab above shows the drift.
Often you do not need samples from p — you need an expectation under p, and you already have samples from somewhere else. The fix is one multiplication per sample.
Suppose the goal is E_p[f(x)] = ∫ f(x)·p(x) dx, but your samples come from a different distribution q. Multiply and divide by q inside the integral: ∫ f(x)·(p(x)/q(x))·q(x) dx. That is now an expectation under q of the function f(x)·w(x), with the importance weight w(x) = p(x)/q(x). So you can keep every sample you drew — you just scale its contribution.
Plain English: some neighbourhoods are over-represented in the proposal and some under-represented. The weight corrects for that. If q gives a region ten times too many samples, each of those samples gets one tenth the vote; if q under-samples a region, each sample there votes ten times louder. It is a census reweighting, not a new election.
Draw x from the proposal q — the shifted curve — then estimate E_p[x²] = 1 by weighting each sample by p(x)/q(x). Slide the proposal away from the target and watch the effective sample size collapse.
The weight is w = p(x)/q(x). Good overlap means weights near 1; bad overlap means a handful of samples dominate. PPO clips these weights for exactly this reason.
The identity is a one-line multiply-and-divide, valid as long as q(x) > 0 wherever p(x) > 0. The variance is where the practical advice lives: the estimator averages f·w, so its variance is Var_q(f·w)/N. If q is tiny where the target mass f·p is large, a handful of enormous weights dominate the average — the estimate is still unbiased, but its spread can be enormous.
Note that weights are not probabilities: here they run up to 1.8 at x = 0.9, and they would grow without bound if q under-covered the right tail. That is the failure mode. Self-normalizing divides by the sum of weights, which tames the scale but introduces a small bias; in practice it usually reduces the mean squared error because the variance term shrinks faster.
Watch the effective sample size in the lab: ESS = (Σw)²/Σw² tells you how many equally-weighted samples your weighted set is worth. When the proposal is far from the target, ESS collapses long before N does — a thousand samples can be worth ten.
You drew x ~ q(x) but need E_p[f(x)]. What do you multiply f(x) by?
For a high-dimensional target you can only evaluate up to a constant, build a random walk whose long-run visits are proportional to the target density. That is Markov chain Monte Carlo.
A Markov chain has no memory beyond its current state: where it goes next depends on where it is, not on the whole path that brought it there. A well-designed chain has a stationary distribution — the long-run fraction of time it spends in each region stops changing. MCMC turns that around: choose the moves so the stationary distribution is exactly the target you cannot sample directly.
Metropolis–Hastings is the foundational recipe. From the current state x, propose a new state x′ from a proposal distribution q(x′|x). Compute the acceptance ratio α = [p(x′)·q(x|x′)] / [p(x)·q(x′|x)] and accept with probability min(1, α); otherwise stay put (and record the stay as a sample). Uphill moves are always taken. Downhill moves are taken in proportion to how much probability they give up. The unknown normalizing constant of p appears in both numerator and denominator and cancels.
A random walk proposes Gaussian steps; uphill moves are always taken, downhill moves with probability p(x′)/p(x). Tune the step size and burn-in, then watch the trace and histogram.
Acceptance near 100% means the steps are too small; acceptance near 0 means they are too large. For a 1-D Gaussian target the sweet spot is around 45% acceptance — in high dimensions, about 23%.
The acceptance rule looks arbitrary until you check detailed balance: the probability flow from x to x′ under the stationary distribution p must equal the flow from x′ to x. When both directions match for every pair of states, p is stationary — the chain’s distribution stops changing.
For a standard normal target with a Gaussian proposal of width s, the acceptance probability has a clean closed form: 2·Φ(−s/2). Plug in real widths:
In high dimensions the optimal acceptance rate for a Gaussian proposal settles near 23.4% (Roberts–Gelman–Gilks). That single number is the reason MCMC libraries print acceptance rates at you: acceptance is the only tuning signal you get for free.
Your proposal is symmetric: q(x′|x) = q(x|x′). What is the acceptance ratio?
A language model emits one logit per token in its vocabulary. How those logits become a choice is entirely the sampler’s job — and three knobs do almost all the work.
Softmax converts logits into probabilities: pᵢ = e^(zᵢ/T) / Σⱼ e^(zⱼ/T). The temperature T divides every logit before exponentiating. Temperature sharpens the distribution when T < 1, flattens it when T > 1, and becomes greedy decoding as T → 0. Top-k keeps only the k most probable tokens and renormalizes. Top-p (nucleus sampling) keeps the smallest set of tokens whose probabilities add up to at least p, then renormalizes — so the candidate set shrinks when the model is confident and grows when it is not.
Greedy decoding always takes the argmax and produces the same answer every time. Sampling uniformly over the whole vocabulary produces gibberish. Beam search sits at the other deterministic extreme: it keeps the B most probable sequences, which maximizes likelihood but tends to produce bland, repetitive text — the highest-probability continuation is often the safest one. Sampling methods exist in the middle, and they are what makes generated text feel alive.
Eight candidate tokens. Reshape the distribution, cut the tail, then draw — and watch the sampled counts follow the final bars.
Temperature never changes which token ranks first — only the gaps between odds. Top-k cuts a fixed count; top-p cuts a fixed mass. Greedy decoding is the T → 0 limit of this same picture.
Split the softmax into numerator and denominator: the ratio of two probabilities is independent of the normalizer, so p₁/p₂ = e^(z₁/T)/e^(z₂/T) = e^((z₁ − z₂)/T). Temperature divides every logit gap by T. With a gap of 2: T = 1 gives odds e² ≈ 7.4 : 1; T = 0.5 gives e⁴ ≈ 55 : 1; T = 2 gives e ≈ 2.7 : 1. The ranking of tokens never changes — only the confidence.
The table is the whole story. Lower T concentrates mass on the leaders; top-p adapts the cut-off to the resulting shape: at T = 0.3 the model is so confident that top-p 0.9 keeps a single token, while at T = 1 it keeps four. That is why top-p generally reads better than top-k — k would still allow k − 1 alternatives on a day the model is certain.
A VAE encodes an input into a distribution, samples from it, and decodes. Sampling has no derivative — until you rewrite it so the randomness is a constant input.
A variational autoencoder compresses an input into a latent distribution N(μ, σ²), draws z from it, and decodes z back into data. Training needs gradients to flow from the reconstruction loss into μ and σ. But a direct draw z ~ N(μ, σ²) is not a differentiable function of μ and σ: nudge μ and a different random number comes out, so there is no well-defined slope to follow. Backpropagation stops at the dice roll, and the encoder never learns.
The reparameterization trick: move the randomness to a parameter-free source. Draw ε ~ N(0, 1) once, then set z = μ + σ·ε. The distribution is identical — a shifted, scaled normal is still normal — but now z is an ordinary differentiable function of μ and σ. For a fixed ε, ∂z/∂μ = 1 and ∂z/∂σ = ε, and the chain rule carries the decoder’s gradient back into the encoder. In a real VAE the encoder outputs log(σ²) for numerical stability and exponentiates half of it to get σ.
Why the distribution is unchanged. A linear transform of a normal is normal. Scaling ε by σ gives mean 0 and variance σ²; adding μ shifts the mean to μ and leaves the variance alone. So μ + σ·ε has exactly the distribution of a direct draw.
Why the gradient exists. With ε held fixed, z is ordinary arithmetic. The derivative of the loss L with respect to μ is the loss gradient at z times 1; with respect to σ it is the loss gradient times ε.
Discrete choices need a different trick. For a categorical draw, the same idea is Gumbel-Softmax: add Gumbel noise g = −ln(−ln u) to the log-probabilities and take a softmax with temperature τ instead of a hard argmax. Numeric check with probabilities [0.5, 0.3, 0.2] and uniforms [0.9, 0.5, 0.1]:
Two related ideas round out the toolbox. Stratified sampling splits the space into equal strata and takes one draw per stratum, which can only lower the variance of a Monte Carlo estimate. And every diffusion step is a reparameterized Gaussian sample: x[t−1] = mean + σ_t·z with z ~ N(0, I) — the same move, repeated a thousand times.
Answer before you look. The temperature, top-p and reparameterization questions are exactly the ones that come up when you tune a real generator.
0 / 6 answered · 0 correct
01What does a temperature below 1.0 do to a language model's output distribution?
02What is the key difference between top-k and top-p (nucleus) sampling?
03Why can't you backpropagate through a standard sampling operation z ~ N(mu, sigma²)?
04In Metropolis–Hastings MCMC, what happens if the proposal standard deviation is set much too large?
05In rejection sampling, what happens to the acceptance rate as the dimensionality of the target distribution increases?
06You have samples from q(x) but need an expectation under p(x). What do you do?
Click a card to swap the lazy description for what it actually means.
Four short problems. Try first; a worked answer is one click away.
Solve u = ½ + arctan(x)/π for x: arctan(x) = π(u − ½), so x = tan(π(u − ½)). Numeric checks: u = 0.5 → x = 0 (the median); u = 0.75 → x = tan(π/4) = 1; u = 0.975 → x = tan(0.475π) ≈ 12.7. The sample histogram matches the bell-shaped PDF, but the tails are so heavy that the mean ∫x·f(x)dx diverges: a single u = 0.9999 gives x ≈ 3183, and the running average jumps every time such a draw appears. The median settles at 0; the mean never does. This is why the Cauchy distribution has no expected value, and a good reminder that “it has a PDF” does not imply “it has a mean”.
With q(x) = 1 on [0, 1], the envelope is M = max p(x). The beta density peaks at x = (a−1)/(a+b−2) = 1/5 = 0.2, and B(2,5) = Γ(2)Γ(5)/Γ(7) = (1!·4!)/6! = 24/720 = 1/30. So p(0.2) = 0.2·0.8⁴·30 = 2.4576, and choosing M = 2.46 gives an acceptance rate of 1/M ≈ 40.7%. For 1,000 accepted samples you expect about 1000/0.407 ≈ 2,457 proposals; the algorithm is exact as long as M ≥ 2.4576 everywhere. If you set M = 2.0 the envelope dips below the peak, some proposals from that region are always accepted, and the accepted sample is quietly biased toward x = 0.2.
Write the integral as an expectation: x = πu with u ~ Uniform(0, 1), so I = π·E[sin(πU)]. The estimator is π·mean(sin(πuᵢ)). Its standard deviation is π·√(Var(sin(πU))/N), and Var(sin(πU)) = E[sin²] − (2/π)² = ½ − 4/π² ≈ 0.0947, so sd ≈ 0.967/√N: about 0.031 at N = 1,000, 0.0097 at N = 10,000, 0.0031 at N = 100,000. One run might give 1.972 (error 0.028), 2.008 (0.008), 1.997 (0.003). Each 100× in samples buys 10× in accuracy — a log–log line of slope −½. Individual runs wobble; fitting many runs recovers the slope.
Work with log p = −(x²y² + x² + y² − 8x − 8y)/2; the unknown normalizer cancels in every ratio. Setting both partials to zero gives x(1+y²) = 4 and y(1+x²) = 4. Subtracting yields (x−y)(1−xy) = 0. The branch x = y solves x³ + x − 4 = 0 at x ≈ 1.379 — a saddle, not a mode. The branch xy = 1 gives x + 1/x = 4, i.e. x = 2 ± √3, so the two modes are (3.732, 0.268) and (0.268, 3.732). At either, xy = 1 and x + y = 4, so log p = −(1 + 14 − 32)/2 = 8.5 and the unnormalized density is e^{8.5} ≈ 4915. Run several chains from different starts. Proposal step 0.1 accepts almost everything but crosses the saddle extremely slowly (a random walk of step 0.1 needs on the order of (4.9/0.1)² ≈ 2,400 accepted steps just to travel between modes). Step 3 jumps freely but rejects most proposals. Around 0.5–1.5 the trace visits both modes and the histogram matches the density; always compare chains and discard burn-in.
You do not need to master these here. Each one gets a proper treatment in its own lesson; the one-line meaning is enough to keep reading. Orange dotted underlines in the prose point back to this list.
This lesson is a port of an open course. Everything here traces back to it — and the next step is running the code yourself.
Lesson text adapted from AI Engineering from Scratch (Phase 01, Lesson 16) and the Math Foundations Notebook reference build. The rejection-rain cover, the estimator bias/variance simulator, the inverse-CDF mapper, the importance-weighting bench, the interactive Metropolis–Hastings chain, the next-token sampler, and the reparameterization figure are original to this page. All labs run in your browser; the MCMC demo is a simplified 1-D random-walk model, labeled as such.