Monte Carlo & MCMC

A probabilistic model is useful only through the quantities you can read off it, and almost every such quantity is an integral: an expectation ${\mathbb{E}}_p[f]$, a marginal $\int p(x,h)dh$, a normalizer $Z=\int \tilde{p}(x)dx$. For the high-dimensional distributions deep models induce, these integrals have no closed form and a grid quadrature costs $O(k^d)$, hopeless past a handful of dimensions. Monte Carlo trades the integral for a sample average, and when even sampling the target is intractable, Markov-chain Monte Carlo generates the samples with a chain designed to converge to it.

Monte Carlo estimation

Every quantity below is an expectation of some function $f$ under a density $p$:

The Monte Carlo estimator draws $n$ independent samples $x_1,\dots ,x_n\sim p$ and returns their empirical mean:

Two facts make it the standard tool of probabilistic inference. It is unbiased (its expectation equals the target) and its variance falls as $1/n$, independent of the dimension $d$.

The decisive consequence is the $1/\sqrt{n}$ law: the typical error shrinks like $\sigma /\sqrt{n}$ regardless of $d$. The dimension does not appear in the rate; it enters only through $\sigma$. This is why Monte Carlo avoids the curse of dimensionality that makes grid methods intractable.¹

The central limit theorem upgrades the variance bound to a full sampling distribution: ${\hat{s}}_n\approx \mathcal{N}\left(s,{\sigma }^2/n\right)$ for large $n$, which is what justifies the error bars ${\hat{s}}_n\pm 1.96\hat{\sigma }/\sqrt{n}$ that accompany any honest Monte Carlo estimate.

The geometric picture

The simplest instance is estimating an area (an integral of an indicator). Scatter points uniformly over a bounding box of known area $A$; the fraction that land in the region $R$ estimates $|R|/A$. With $f=1[x\in R]$ and $p$ uniform on the box, ${\hat{s}}_n$ reduces to that landing fraction: Monte Carlo integration with $f$ an indicator.

Importance sampling

Often we cannot draw from $p$, or $f$ concentrates its mass where $p$ rarely puts samples. To address this, sample from a tractable proposal $q$ and correct the mismatch by reweighting. Multiply and divide by $q$:

Sampling $x_i\sim q$ and averaging gives the importance-sampling estimator with weights $w_i=p(x_i)/q(x_i)$:

The estimator is unbiased for any $q$ with $q(x)>0$ wherever $p(x)f(x)\ne 0$, but its variance depends sharply on the choice of $q$. Minimizing ${\mathrm{Var}}_q[wf]$ over $q$ yields the optimal proposal.

The optimum is unreachable in practice, since it presupposes the very integral $s$ we are after (it needs $\int p|f|$). It still dictates the design rule: put proposal mass where $p\left|f\right|$ is large. The corresponding failure is a proposal with too little mass in those regions.

Markov chains

When the target $p(x)$ admits no direct sampler, the typical case for the energy-based models of the next chapter, whose density is $p(x)=\tilde{p}(x)/Z$ with an intractable $Z$, we abandon iid draws and build a Markov chain whose long-run distribution is $p$. A chain is governed by a transition kernel $T(x^{\prime }\mid x)$, the probability of moving to $x^{\prime }$ given the current state $x$. Starting from $q^{(0)}$, the distribution evolves by one matrix–vector product per step:

For example, take three states $A,B,C$ with the transition matrix below, where entry $T_{ij}$ is the probability of moving from state $j$ to state $i$ (columns sum to $1$). One step of the chain is the matrix–vector product $q^{(t+1)}=T{q}^{(t)}$; iterating it drives any starting distribution toward the stationary $p$, the eigenvector of $T$ with eigenvalue $1$.

The transition matrix and its fixed point are

The chain forgets its start because every column mixes probability across states; only the eigenvector direction persists under repeated multiplication. MCMC constructs such a $T$ so that its fixed point $p$ is a distribution chosen in advance.

A convenient sufficient condition for stationarity is detailed balance, which we will use to design kernels with a prescribed stationary distribution.

Stationarity alone does not guarantee the chain reaches $p$ from an arbitrary start. That requires ergodicity: the chain must be irreducible (every state reachable from every other) and aperiodic. Under ergodicity the stationary distribution is unique and $q^{(t)}\to p$ from any $q^{(0)}$, at a rate set by the second-largest eigenvalue $\left|{\lambda }_2\right|$ of $T$.³

Markov-chain Monte Carlo

MCMC designs a kernel $T$ whose stationary distribution equals the target $p$, runs the chain, and treats the visited states as (correlated) samples. Two constructions dominate.

Metropolis–Hastings

Pick any proposal $g(x^{\prime }\mid x)$ and make every proposed move conditional on an accept/reject test built from the density ratio. The key property is that the test needs $p$ only up to its normalizer, since $Z$ cancels in the ratio.

When the proposal is symmetric, $g(x^{\prime }\mid x)=g(x\mid x^{\prime })$ (e.g. a Gaussian random walk), the $g$ factors cancel and the rule collapses to the original Metropolis acceptance $A=\min (1,\tilde{p}(x^{\prime })/\tilde{p}(x))$: always accept a move uphill in density, accept a downhill move with probability equal to the density ratio.

For example, suppose a symmetric random walk sits at $x$ with $\tilde{p}(x)=2.0$ and proposes $x^{\prime }$ with $\tilde{p}(x^{\prime })=0.8$. The ratio is $0.8/2.0=0.4$, so $A=\min (1,0.4)=0.4$: draw $u\sim \text{Uniform}(0,1)$ and take the move only if $u\le 0.4$, otherwise stay at $x$ and record it again. Had the proposal climbed to $\tilde{p}(x^{\prime })=3.0$, the ratio $1.5$ would clip to $A=1$ and the move would always be taken. The chain therefore spends time in each region in proportion to $\tilde{p}$, satisfying the stationary condition $Z$ cancelled out of.

Gibbs sampling

When the full conditionals $p(x_i\mid x_{-i})$, the distribution of one coordinate given all the others, are tractable, Gibbs sampling dispenses with acceptance entirely: it sweeps the coordinates, resampling each from its conditional. This is the natural sampler for energy-based models, where each unit's conditional is a simple sigmoid of its neighbors.⁴

Each single-coordinate update leaves $p$ invariant, and a full sweep (a composition of invariant kernels) therefore does too.

The reason Gibbs sidesteps the acceptance test entirely is that the full conditional, like the M–H ratio, never sees $Z$. Write the joint as $p(x)=\tilde{p}(x)/Z$ and condition on the fixed coordinates $x_{-i}$:

The global normalizer $Z$ cancels top and bottom, leaving only a sum over the single coordinate $x_i$ — a one-dimensional normalization that is tractable even when the $d$-dimensional $Z$ is not. For a binary energy-based unit this collapses to a sigmoid: if $\tilde{p}(x_i,x_{-i})=\exp \left(x_i{b}_i+x_i{\sum }_j{W}_{ij}{x}_j\right)$ with $x_i\in \{0,1\}$, then $p(x_i=1\mid x_{-i})=\sigma \left(b_i+{\sum }_j{W}_{ij}{x}_j\right)$, which is why Gibbs is the native sampler of Boltzmann machines.

Geometrically, Gibbs moves are axis-aligned: each step slides along one coordinate to a value drawn from the conditional slice, tracing a staircase through the density.

Mixing

MCMC samples are correlated (each state is a perturbation of the last), so the practical question is how many steps separate two effectively independent draws. Two diagnostics govern this.

The variance penalty is the whole cost of correlation. Expanding the variance of the chain average $\hat{s}=\frac{1}{N}{\sum }_{t}f(x^{(t)})$ keeps the cross terms that iid sampling would have killed:

where the last step takes $N\gg \tau$ so the edge correction $1-k/N$ drops out and the sum telescopes into the integrated autocorrelation time $\tau =1+2{\sum }_{k\ge 1}\rho (k)$. Since $\rho (k)\approx {\left|{\lambda }_2\right|}^k$ is geometric, $\tau \approx (1+\left|{\lambda }_2\right|)/(1-\left|{\lambda }_2\right|)$, which blows up as $\left|{\lambda }_2\right|\to 1$. Thus $N$ MCMC steps are worth only $N/\tau$ iid samples. A chain mixes quickly when $\left|{\lambda }_2\right|$ is small ($\tau$ small) and slowly when $\left|{\lambda }_2\right|\to 1$.

The pathological case is a multimodal target. To travel between two modes of $p(x)=\tilde{p}(x)/Z$, a local sampler must cross a region where $\tilde{p}$, and hence the acceptance probability, is exponentially small. The chain is trapped in one mode for an exponentially long time: the energy barrier between modes is the geometric cause of slow mixing.

Heavily correlated coordinates slow Gibbs for the same reason: each axis-aligned move is confined to a thin conditional slice, so the staircase advances in small steps along the ridge of the density.

The partition function

The motivating problem of Goodfellow Ch. 18 is learning an energy-based model $p(x)=\tilde{p}(x)/Z$, where $\tilde{p}(x)=\exp (-E(x))$ is tractable but the normalizer $Z={\sum }_x\tilde{p}(x)$ is an intractable sum over all configurations. Maximum-likelihood training needs the gradient of $\log p$, and that gradient splits into two terms: one cheap, one requiring samples from the model itself.

The second term is itself an expectation under the model, which is the crux:

The two phases have a physical reading: the positive phase pushes down the energy of observed data, while the negative phase pushes up the energy of the model's own fantasies, sampled by a Markov chain. Learning halts when the two balance, when the model's samples are statistically indistinguishable from the data.⁵

Annealed importance sampling

Detecting convergence or comparing two trained models still needs $Z$ itself: for a held-out likelihood you must normalize. A naive importance-sampling estimate $Z\approx \frac{1}{n}\sum \tilde{p}(x_i)/q(x_i)$ collapses because a tractable $q$ never overlaps a sharp model. Annealed importance sampling (AIS) addresses this by bridging from a simple base $p_0$ (with known $Z_0$) to the target $p_T=p$ through a sequence of intermediate distributions $p_0,p_1,\dots ,p_T$, typically the geometric path $p_{\beta }(x)\propto p_0(x{)}^{1-\beta }\tilde{p}(x{)}^{\beta }$ for $0={\beta }_0<\cdots <{\beta }_T=1$.

Each rung's two distributions overlap well, so each ratio is low-variance; the anneal schedule converts one intractable importance-sampling problem into many easy ones, the same divide-and-conquer logic that lets simulated annealing cross the energy barriers that trap a single-temperature chain.⁶

Gradient-guided samplers and the road to diffusion

The random-walk Metropolis of this lesson proposes moves blindly, so in high dimensions its acceptance rate forces tiny steps and the chain moves slowly. The samplers that dominate modern practice fix this by proposing moves along the gradient of the log-density — the same ${\nabla }_x\log \tilde{p}(x)$ that appears nowhere in Metropolis but everywhere in deep learning.

Hamiltonian Monte Carlo (HMC; Neal, 2011, in the Handbook of MCMC) augments the state $x$ with a momentum $r$ and simulates physical dynamics under the potential energy $E(x)=-\log \tilde{p}(x)$: a particle given a random kick coasts along the density's contours for many steps before a single Metropolis accept/reject corrects the numerical error. Because the proposal follows the gradient rather than diffusing, HMC takes long, high-acceptance moves and its autocorrelation time barely grows with dimension. NUTS (the No-U-Turn Sampler; Hoffman & Gelman, 2014, JMLR) removes HMC's one awkward tuning knob — the trajectory length — by running the simulation until it starts to double back, and is the default engine inside probabilistic-programming systems like Stan and PyMC.

Langevin dynamics is the one-step limit of the same idea: propose $x^{\prime }=x+\frac{{ϵ}^2}{2}{\nabla }_x\log \tilde{p}(x)+ϵz$ with $z\sim \mathcal{N}(0,I)$ — a gradient-ascent step on the log-density plus calibrated noise. This is the bridge to the generative models of the next module. Score-based generative models (Song & Ermon, 2019, NeurIPS) and the equivalent denoising diffusion models (Ho et al., 2020, NeurIPS) learn the score ${\nabla }_x\log p(x)$ with a neural network and then sample by running exactly this Langevin update, annealed across noise levels in the same spirit as AIS. The negative-phase problem that made energy-based learning hard — needing samples from the model — is answered by learning the score directly and letting gradient-guided Langevin sampling do the rest, which is why diffusion models train stably where Boltzmann machines did not.

Takeaways

• Monte Carlo estimates any expectation $s={\mathbb{E}}_p[f]$ by the sample mean ${\hat{s}}_n=\frac{1}{n}{\sum }_{i}f(x_i)$: unbiased, with variance ${\sigma }^2/n$ and error $\sigma /\sqrt{n}$, a rate independent of dimension, which is why it beats grid quadrature in high $d$.
• Importance sampling draws from a tractable proposal $q$ and reweights by $w=p/q$; the optimal proposal is $q^{\star }\propto p\left|f\right|$, and a $q$ with lighter tails than $p$ blows the variance up.
• A Markov chain with kernel $T$ converges to a unique stationary $p$ when it is ergodic; detailed balance is the design tool that fixes $p$ as the stationary distribution.
• MCMC builds such a chain: Metropolis–Hastings accepts a proposed move with probability $\min (1,\tilde{p}(x^{\prime })g(x\mid x^{\prime })/(\tilde{p}(x)g(x^{\prime }\mid x)))$ — $Z$ cancels; Gibbs resamples each coordinate from its full conditional, leaving $p$ invariant.
• Mixing is governed by $\left|{\lambda }_2\right|$: burn-in discards the transient, autocorrelation time $\tau$ shrinks the effective sample size to $N/\tau$, and multimodal targets mix slowly because energy barriers trap local chains.
• The partition function $Z$ makes energy-based learning hard: the log-likelihood gradient is a positive phase (data) minus a negative phase (model samples via MCMC); annealed importance sampling estimates $Z$ itself by bridging from a tractable base. See energy-based and Boltzmann machines.