---
title: Monte Carlo & MCMC
module: Probabilistic Methods
moduleNumber: 8
lessonNumber: 2
order: 802
summary: >
  Most quantities of interest in a probabilistic model are integrals nobody can
  compute in closed form: expectations, marginals, partition functions. Monte
  Carlo replaces the integral with an average over samples; importance sampling
  reweights samples from a tractable proposal; and when even sampling the target
  is hard, Markov-chain Monte Carlo builds a chain whose stationary distribution
  _is_ the target. We derive Metropolis–Hastings and Gibbs, analyze mixing, and
  close on the partition-function gradient that powers energy-based learning.
topics: [Probabilistic Methods]
sources:
  - book: Goodfellow
    ref: "Ch. 17 — Monte Carlo Methods"
  - book: Goodfellow
    ref: "Ch. 18 — Confronting the Partition Function"
---

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$:

$$
s = \mathbb{E}_{x \sim p}[f(x)] = \int p(x)\,f(x)\,dx
\quad\text{or}\quad
\sum_x p(x)\,f(x).
$$

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

$$
\hat s_n = \frac{1}{n}\sum_{i=1}^{n} f(x_i),
\qquad x_i \overset{\text{iid}}{\sim} p.
$$

> **Definition (Monte Carlo estimator).** Given iid draws $x_i \sim p$, the
> estimator $\hat s_n = \tfrac1n\sum_i f(x_i)$ of $s = \mathbb{E}_p[f]$. It
> converts the analytic problem of integration into the statistical problem of
> averaging samples.

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$.

> **Theorem (Unbiasedness and variance).** If $x_i \overset{\text{iid}}{\sim} p$
> and $\sigma^2 = \Var_p[f(x)] < \infty$, then
> $\mathbb{E}[\hat s_n] = s$ and $\Var[\hat s_n] = \sigma^2 / n$.

> **Proof.** By linearity, $\mathbb{E}[\hat s_n] = \tfrac1n\sum_i \mathbb{E}[f(x_i)]
> = \tfrac1n \cdot n\,s = s$, so the estimator is unbiased. For the variance,
> independence kills all cross terms, leaving
> $$
> \Var[\hat s_n]
> = \frac{1}{n^2}\sum_{i=1}^{n}\Var[f(x_i)]
> = \frac{1}{n^2}\cdot n\,\sigma^2
> = \frac{\sigma^2}{n}.
> $$
> The standard error is therefore $\sigma/\sqrt n$. $\qed$

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.[^gf-mc]

| Method | Cost for error $\varepsilon$ | Scaling in dimension $d$ |
| --- | --- | --- |
| Grid / deterministic quadrature | $O(\varepsilon^{-d})$ | exponential in $d$ |
| Monte Carlo | $O(\sigma^2 \varepsilon^{-2})$ | _no explicit $d$_ — only via $\sigma^2$ |

The central limit theorem upgrades the variance bound to a full sampling
distribution: $\hat s_n \approx \mathcal{N}\!\parens{s, \sigma^2/n}$ 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 = \mathbf{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.

$$
% caption: Monte Carlo area estimation: the fraction of uniform points landing inside region $R$ (green hits) estimates its area, with error shrinking as $1/\sqrt n$.
\begin{tikzpicture}[>=stealth, font=\small, scale=2.4]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % bounding box
  \draw[black, thick] (0,0) rectangle (1,1);
  \node[font=\footnotesize, black, anchor=south east] at (1,1.03) {box, area $1$};
  % the curve y = x^2 bounding region R below it... use quarter circle region instead
  % region R = under curve y = sqrt(1 - (1-x)^2)? keep simple: under y = x curve? use quarter disk
  \draw[acc, very thick] plot[domain=0:1, samples=40] (\x, {sqrt(1-(1-\x)*(1-\x))});
  \node[acc, font=\footnotesize, anchor=north] at (0.5,-0.03) {region $R$};
  % HIT points (inside quarter disk centered at (1,0) radius 1): y < sqrt(1-(1-x)^2)
  \foreach \px/\py in {0.20/0.10, 0.45/0.25, 0.30/0.40, 0.62/0.30, 0.55/0.55, 0.78/0.20, 0.40/0.12, 0.70/0.50, 0.85/0.35, 0.50/0.08}
    \fill[green] (\px,\py) circle (0.014);
  % MISS points (outside): above the curve
  \foreach \px/\py in {0.15/0.80, 0.35/0.90, 0.10/0.55, 0.25/0.70, 0.50/0.85, 0.40/0.75, 0.20/0.92, 0.30/0.62, 0.08/0.78, 0.18/0.68}
    \fill[black] (\px,\py) circle (0.014);
  \node[green, font=\footnotesize, anchor=west] at (1.08,0.62) {hit ($f=1$)};
  \node[black, font=\footnotesize, anchor=west] at (1.08,0.42) {miss ($f=0$)};
  \node[font=\footnotesize, anchor=west] at (1.08,0.22) {$\hat s_n=\dfrac{\#\text{hits}}{n}$};
\end{tikzpicture}
$$

## 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$:

$$
s = \mathbb{E}_{x \sim p}[f(x)]
= \int p(x)\,f(x)\,dx
= \int q(x)\,\frac{p(x)}{q(x)}\,f(x)\,dx
= \mathbb{E}_{x \sim q}\!\brackets{\frac{p(x)}{q(x)}\,f(x)}.
$$

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

$$
\hat s_q = \frac{1}{n}\sum_{i=1}^{n} w_i\,f(x_i),
\qquad w_i = \frac{p(x_i)}{q(x_i)}.
$$

> **Definition (Importance weight).** The ratio $w(x) = p(x)/q(x)$ correcting a
> sample drawn from the proposal $q$ back to the target $p$. When only the
> unnormalized $\tilde p = Z\,p$ is available, the **self-normalized** estimator
> uses $\tilde w_i = \tilde p(x_i)/q(x_i)$ and divides by $\sum_j \tilde w_j$,
> sidestepping $Z$.

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
$\Var_q[w\,f]$ over $q$ yields the optimal proposal.

> **Theorem (Optimal proposal).** The variance of $\hat s_q$ is minimized by
> $q^\star(x) \propto p(x)\,\abs{f(x)}$. When $f \ge 0$, this proposal gives
> _zero_ variance.

> **Proof.** Write $g(x) = p(x)f(x)/q(x)$, so $\hat s_q$ averages $g(x_i)$ with
> $x_i \sim q$ and $\mathbb{E}_q[g] = s$. Then
> $\Var_q[g] = \mathbb{E}_q[g^2] - s^2$, and only the first term
> depends on $q$:
> $$
> \mathbb{E}_q[g^2] = \int q(x)\,\frac{p(x)^2 f(x)^2}{q(x)^2}\,dx
> = \int \frac{p(x)^2 f(x)^2}{q(x)}\,dx.
> $$
> By the Cauchy–Schwarz inequality applied to
> $\int \tfrac{p|f|}{\sqrt q}\cdot\sqrt q\,dx$,
> $$
> \parens{\int p(x)\,\abs{f(x)}\,dx}^2
> \le \parens{\int \frac{p^2 f^2}{q}\,dx}\parens{\int q\,dx}
> = \mathbb{E}_q[g^2],
> $$
> with equality iff $q(x) \propto p(x)\abs{f(x)}$. That choice attains the
> lower bound, and for $f \ge 0$ the bound equals $s^2$, so
> $\Var = s^2 - s^2 = 0$. $\qed$

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\,\abs{f}$ is large.** The corresponding failure is a proposal with too
little mass in those regions.

> **Remark (Variance blow-up).** If $q$ has lighter tails than $p$, the ratio
> $w = p/q \to \infty$ in regions $q$ under-samples. A rare draw there carries an
> enormous weight, $\Var_q[w] = \infty$ is possible, and a handful of
> samples dominate the average: the estimator becomes high-variance and unreliable
> despite remaining unbiased.[^gf-importance]

$$
% caption: Importance sampling: samples drawn from proposal $q$ (black) are reweighted by $w=p/q$ to estimate under the target $p$ (blue).
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % axes
  \draw[->, thick] (-3.4,0) -- (3.7,0) node[right, font=\footnotesize] {$x$};
  \draw[->, thick] (0,-0.2) -- (0,2.7) node[above, font=\footnotesize] {density};
  % target p: gaussian shifted right, peaked
  \draw[acc, very thick] plot[domain=-3.2:3.5, samples=90]
    (\x, {2.3*exp(-(\x-0.9)*(\x-0.9)/0.5)});
  \node[acc, font=\footnotesize, anchor=west] at (1.55,2.05) {target $p$};
  % proposal q: broader gaussian centered left
  \draw[black, very thick] plot[domain=-3.2:3.5, samples=90]
    (\x, {1.5*exp(-(\x+0.2)*(\x+0.2)/1.6)});
  \node[black, font=\footnotesize, anchor=east] at (-1.5,1.35) {\texttt{proposal} $q$};
  % weight annotations
  \draw[->, green, thick] (1.6,0.55) -- (1.6,1.55);
  \node[green, font=\footnotesize, anchor=west] at (1.66,1.1) {$w>1$};
  \draw[->, red, thick] (-1.4,0.95) -- (-1.4,0.35);
  \node[red, font=\footnotesize, anchor=east] at (-1.46,0.7) {$w<1$};
\end{tikzpicture}
$$

## Markov chains

When the target $p(x)$ admits no direct sampler, the typical case for the
[energy-based models](/deep-learning/generative-models/energy-based-and-boltzmann-machines)
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' \mid x)$, the
probability of moving to $x'$ given the current state $x$. Starting from $q^{(0)}$,
the distribution evolves by one matrix–vector product per step:

$$
q^{(t+1)}(x') = \sum_{x} T(x' \mid x)\,q^{(t)}(x),
\qquad
q^{(t)} = T^t\,q^{(0)}.
$$

> **Definition (Stationary distribution).** A distribution $p$ is **stationary**
> for kernel $T$ if it is a fixed point of the update,
> $p(x') = \sum_x T(x' \mid x)\,p(x)$ — once the chain reaches $p$, it stays there.
> In matrix terms $p$ is the eigenvector of $T$ with eigenvalue $1$.

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$.

$$
% caption: A three-state Markov chain. Each arrow $j \to i$ carries the transition probability $T(i \mid j)$; self-loops hold the state. The stationary $p$ satisfies $p = Tp$.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \tikzset{st/.style={circle, draw=black, thick, minimum size=1.0cm, font=\normalsize}}
  \node[st] (A) at (0,0) {A};
  \node[st] (B) at (4,0) {B};
  \node[st] (C) at (2,-2.6) {C};
  % edges with probabilities
  \draw[->, acc, thick] (A) to[bend left=18] node[above, font=\scriptsize, black]{0.6} (B);
  \draw[->, acc, thick] (B) to[bend left=18] node[below, font=\scriptsize, black]{0.3} (A);
  \draw[->, acc, thick] (B) to[bend left=18] node[right, font=\scriptsize, black]{0.5} (C);
  \draw[->, acc, thick] (C) to[bend left=18] node[left, font=\scriptsize, black]{0.4} (B);
  \draw[->, acc, thick] (A) to[bend right=18] node[left, font=\scriptsize, black]{0.2} (C);
  \draw[->, acc, thick] (C) to[bend right=18] node[right, font=\scriptsize, black]{0.3} (A);
  % self-loops
  \draw[->, acc, thick] (A) to[out=150, in=210, looseness=6] node[left, font=\scriptsize, black]{0.2} (A);
  \draw[->, acc, thick] (B) to[out=30, in=-30, looseness=6] node[right, font=\scriptsize, black]{0.2} (B);
  \draw[->, acc, thick] (C) to[out=250, in=290, looseness=6] node[below, font=\scriptsize, black]{0.3} (C);
\end{tikzpicture}
$$

The transition matrix and its fixed point are

$$
T =
\begin{pmatrix}
0.2 & 0.3 & 0.3\\
0.6 & 0.2 & 0.4\\
0.2 & 0.5 & 0.3
\end{pmatrix},
\qquad
p = Tp,
\qquad
p \approx (0.27,\ 0.40,\ 0.33).
$$

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.

> **Theorem (Detailed balance $\Rightarrow$ stationarity).** If $T$ and $p$ satisfy
> $p(x)\,T(x' \mid x) = p(x')\,T(x \mid x')$ for all $x, x'$, then $p$ is stationary
> for $T$.

> **Proof.** Sum the balance equation over $x$:
> $$
> \sum_x p(x)\,T(x' \mid x)
> = \sum_x p(x')\,T(x \mid x')
> = p(x')\sum_x T(x \mid x')
> = p(x'),
> $$
> using $\sum_x T(x \mid x') = 1$ (the kernel is a valid distribution over the next
> state). The left side is the one-step update of $p$, which equals $p(x')$ — so $p$
> is a fixed point. $\qed$

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 $\abs{\lambda_2}$ of $T$.[^gf-markov]

| Property | Meaning | Why it matters |
| --- | --- | --- |
| Stationary | $p$ is a fixed point of $T$ | the chain leaves $p$ invariant |
| Irreducible | all states communicate | no region is trapped or unreachable |
| Aperiodic | no forced cyclic return | the chain mixes rather than oscillates |
| Ergodic | irreducible $+$ aperiodic | $q^{(t)} \to p$ uniquely, from any start |

$$
% caption: A Markov chain converging to its stationary distribution $p$ (dashed): the state histogram (blue) approaches $p$ as the burn-in transient decays at rate $\abs{\lambda_2}^t$.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % three mini-panels: t=0, t=10, t=200
  \foreach \px/\lab/\h in {0/{$t=0$}/0, 4.4/{$t=10$}/1, 8.8/{$t=200$}/2} {
    \begin{scope}[xshift=\px cm]
      \draw[->, thick] (-0.1,0) -- (3.3,0) node[right, font=\scriptsize] {$x$};
      \draw[->, thick] (0,-0.1) -- (0,2.4);
      \node[font=\footnotesize, anchor=south] at (1.5,2.3) {\lab};
    \end{scope}
  }
  % target curve (same on all three): bell over [0,3] peak ~1.8 at x=1.5
  \foreach \px in {0,4.4,8.8} {
    \begin{scope}[xshift=\px cm]
      \draw[black, dashed, thick] plot[domain=0.1:3.1, samples=40]
        (\x, {1.8*exp(-(\x-1.5)*(\x-1.5)/0.6)});
    \end{scope}
  }
  % t=0 histogram: all mass on the left (far from target)
  \begin{scope}[xshift=0cm]
    \draw[acc, very thick, fill=acc!15] (0.15,0) rectangle (0.55,2.0);
    \draw[acc, very thick, fill=acc!15] (0.55,0) rectangle (0.95,1.1);
    \draw[acc, very thick, fill=acc!15] (0.95,0) rectangle (1.35,0.4);
  \end{scope}
  % t=10 histogram: spreading toward target, still skewed
  \begin{scope}[xshift=4.4cm]
    \draw[acc, very thick, fill=acc!15] (0.35,0) rectangle (0.75,0.9);
    \draw[acc, very thick, fill=acc!15] (0.75,0) rectangle (1.15,1.6);
    \draw[acc, very thick, fill=acc!15] (1.15,0) rectangle (1.55,1.7);
    \draw[acc, very thick, fill=acc!15] (1.55,0) rectangle (1.95,1.2);
    \draw[acc, very thick, fill=acc!15] (1.95,0) rectangle (2.35,0.6);
  \end{scope}
  % t=200 histogram: matches target bell
  \begin{scope}[xshift=8.8cm]
    \draw[acc, very thick, fill=acc!15] (0.35,0) rectangle (0.75,0.5);
    \draw[acc, very thick, fill=acc!15] (0.75,0) rectangle (1.15,1.25);
    \draw[acc, very thick, fill=acc!15] (1.15,0) rectangle (1.55,1.75);
    \draw[acc, very thick, fill=acc!15] (1.55,0) rectangle (1.95,1.7);
    \draw[acc, very thick, fill=acc!15] (1.95,0) rectangle (2.35,1.15);
    \draw[acc, very thick, fill=acc!15] (2.35,0) rectangle (2.75,0.45);
  \end{scope}
  \node[black, font=\footnotesize, anchor=west] at (8.95,2.05) {target $p$};
\end{tikzpicture}
$$

## 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' \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.

$$
A(x' \mid x)
= \min\!\parens{1,\;
\frac{p(x')\,g(x \mid x')}{p(x)\,g(x' \mid x)}}
= \min\!\parens{1,\;
\frac{\tilde p(x')\,g(x \mid x')}{\tilde p(x)\,g(x' \mid x)}}.
$$

```algorithm
caption: $\textsc{MetropolisHastings}(\tilde p, g, x_0, N)$ — sample $p \propto \tilde p$
initialize $x \gets x_0$
for $t \gets 1$ to $N$ do
  propose $x' \sim g(\cdot \mid x)$
  $a \gets \min\!\parens{1,\ \tilde p(x')\,g(x \mid x') \,/\, (\tilde p(x)\,g(x' \mid x))}$ // acceptance ratio
  draw $u \sim \text{Uniform}(0,1)$
  if $u \le a$ then // accept
    $x \gets x'$
  else // reject, stay put
    $x \gets x$
  record $x$ as sample $x^{(t)}$
return $x^{(1)}, \dots, x^{(N)}$
```

> **Theorem (M–H invariance).** The Metropolis–Hastings kernel
> $T(x' \mid x) = g(x' \mid x)\,A(x' \mid x)$ (for $x' \ne x$) satisfies detailed
> balance with $p$, hence leaves $p$ invariant.

> **Proof.** Take $x \ne x'$ and assume WLOG the acceptance ratio is $\le 1$ for the
> forward move, so $A(x' \mid x) = \tfrac{p(x')g(x \mid x')}{p(x)g(x' \mid x)}$ and
> the reverse $A(x \mid x') = 1$. Then
> $$
> p(x)\,T(x' \mid x)
> = p(x)\,g(x' \mid x)\,\frac{p(x')\,g(x \mid x')}{p(x)\,g(x' \mid x)}
> = p(x')\,g(x \mid x'),
> $$
> while $p(x')\,T(x \mid x') = p(x')\,g(x \mid x')\cdot 1$. The two sides match, so
> detailed balance holds; by the previous theorem $p$ is stationary. $\qed$

When the proposal is symmetric, $g(x' \mid x) = g(x \mid x')$ (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')/\tilde p(x))$: always accept a move uphill in
density, accept a downhill move with probability equal to the density ratio.

$$
% caption: Metropolis on a target $\tilde p$ with a symmetric proposal. An uphill proposal ($\tilde p(x') > \tilde p(x)$) is always accepted; a downhill one is accepted with probability $\tilde p(x')/\tilde p(x)$, so the chain still explores the low-density valleys.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  \draw[->, thick] (-0.2,0) -- (7.0,0) node[right, font=\footnotesize] {$x$};
  \draw[->, thick] (0,-0.2) -- (0,3.0);
  % target density (single broad hump)
  \draw[acc, very thick] plot[domain=0.2:6.6, samples=120]
    (\x, {2.4*exp(-(\x-3.3)*(\x-3.3)/3.2)});
  \node[acc, font=\footnotesize, anchor=west] at (4.85,2.5) {target $\tilde p$};
  % current state x
  \fill[black] (2.2,1.72) circle (2.3pt);
  \node[black, font=\footnotesize, anchor=south east] at (2.2,1.78) {$x$};
  \draw[black, dashed] (2.2,0) -- (2.2,1.72);
  % uphill proposal (accepted)
  \fill[green] (3.1,2.32) circle (2.3pt);
  \draw[green, dashed] (3.1,0) -- (3.1,2.32);
  \node[green, font=\scriptsize, anchor=south] at (3.1,2.4) {uphill: accept};
  \draw[->, green, thick] (2.35,1.72) to[bend left=12] (2.95,2.28);
  % downhill proposal (maybe)
  \fill[red] (5.4,1.05) circle (2.3pt);
  \draw[red, dashed] (5.4,0) -- (5.4,1.05);
  \node[red, font=\scriptsize, anchor=west] at (5.5,1.1) {downhill:};
  \node[red, font=\scriptsize, anchor=west] at (5.5,0.72) {accept w.p. ratio};
  \draw[->, red, thick] (2.35,1.68) to[bend right=10] (5.25,1.1);
\end{tikzpicture}
$$

For example, suppose a symmetric random walk
sits at $x$ with $\tilde p(x) = 2.0$ and proposes $x'$ with $\tilde p(x') = 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') = 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.[^gf-gibbs]

```algorithm
caption: $\textsc{Gibbs}(p, x_0, N)$ — coordinate-wise resampling from full conditionals
initialize $x \gets x_0 = (x_1, \dots, x_d)$
for $t \gets 1$ to $N$ do
  for $i \gets 1$ to $d$ do // one sweep over coordinates
    $x_i \gets$ sample from $p(x_i \mid x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d)$
  record $x$ as sample $x^{(t)}$
return $x^{(1)}, \dots, x^{(N)}$
```

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

> **Theorem (Gibbs invariance).** Resampling coordinate $i$ from its full conditional
> $p(x_i \mid x_{-i})$ leaves the joint $p$ invariant.

> **Proof.** The update kernel for the $i$-th step is
> $T_i(x' \mid x) = p(x_i' \mid x_{-i})\,\mathbf{1}[x_{-i}' = x_{-i}]$: it changes
> only $x_i$ and draws the new value from the conditional. Apply it to a state
> distributed as $p$ and marginalize the old $x_i$:
> $$
> \sum_{x_i} p(x_i, x_{-i})\,p(x_i' \mid x_{-i})
> = p(x_i' \mid x_{-i}) \sum_{x_i} p(x_i, x_{-i})
> = p(x_i' \mid x_{-i})\,p(x_{-i}).
> $$
> The factorization $p(x_i' \mid x_{-i})\,p(x_{-i}) = p(x_i', x_{-i})$ is the chain
> rule, so the result is $p$ evaluated at the new state — the kernel reproduces $p$.
> A sweep $T = T_d \circ \cdots \circ T_1$ is a composition of $p$-invariant kernels,
> hence $p$-invariant. $\qed$

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}$:

$$
p(x_i \mid x_{-i})
= \frac{p(x_i, x_{-i})}{p(x_{-i})}
= \frac{\tilde p(x_i, x_{-i})/Z}{\sum_{x_i'} \tilde p(x_i', x_{-i})/Z}
= \frac{\tilde p(x_i, x_{-i})}{\sum_{x_i'} \tilde p(x_i', 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\!\parens{x_i\,b_i + x_i \sum_{j} W_{ij} x_j}$ with
$x_i \in \{0,1\}$, then $p(x_i = 1 \mid x_{-i}) = \sigma\!\parens{b_i + \sum_j W_{ij} x_j}$,
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.

$$
% caption: Gibbs sampling on a 2D density resamples one coordinate at a time, giving an axis-aligned staircase path; strong correlation forces tiny steps and slow mixing.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  % axes
  \draw[->, thick] (-0.2,0) -- (5.2,0) node[right, font=\footnotesize] {$x_1$};
  \draw[->, thick] (0,-0.2) -- (0,4.6) node[above, font=\footnotesize] {$x_2$};
  % correlated density: tilted concentric ellipses
  \foreach \r in {0.5,1.0,1.5} {
    \draw[black, rotate around={38:(2.6,2.3)}] (2.6,2.3) ellipse ({\r*2.0} and \r);
  }
  \node[black, font=\footnotesize, anchor=west] at (3.7,4.2) {density};
  % staircase path: alternate horizontal (resample x1) and vertical (resample x2)
  \coordinate (a) at (1.1,1.2);
  \coordinate (b) at (2.5,1.2);
  \coordinate (c) at (2.5,2.0);
  \coordinate (d) at (3.4,2.0);
  \coordinate (e) at (3.4,2.7);
  \coordinate (f) at (4.0,2.7);
  \draw[acc, very thick, ->] (a) -- (b);
  \draw[acc, very thick, ->] (b) -- (c);
  \draw[acc, very thick, ->] (c) -- (d);
  \draw[acc, very thick, ->] (d) -- (e);
  \draw[acc, very thick, ->] (e) -- (f);
  \node[acc, font=\scriptsize, anchor=north] at (1.8,0.75) {resample $x_1$};
  \node[acc, font=\scriptsize, anchor=west] at (4.05,1.35) {resample $x_2$};
  \fill[green] (a) circle (2.2pt);
  \node[green, font=\footnotesize, anchor=east] at (1.0,1.2) {start};
  \fill[acc] (f) circle (2.2pt);
\end{tikzpicture}
$$

## 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.

> **Definition (Burn-in).** The initial transient before the chain reaches its
> stationary regime. The early samples reflect the arbitrary start $q^{(0)}$, not $p$,
> and are discarded. Its length is set by the relaxation time
> $\tau_{\text{relax}} \approx 1/(1 - \abs{\lambda_2})$.

> **Definition (Autocorrelation).** The correlation $\rho(k) =
> \Corr(f(x^{(t)}), f(x^{(t+k)}))$ between samples $k$ steps apart.
> It decays geometrically, $\rho(k) \approx \abs{\lambda_2}^{k}$; the
> **integrated autocorrelation time** $\tau = 1 + 2\sum_{k\ge 1}\rho(k)$ inflates
> the variance of the chain estimate to $\Var[\hat s] \approx
> \tau\,\sigma^2/N$, so the **effective sample size** is $N_{\text{eff}} = N/\tau$.

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

$$
\Var[\hat s]
= \frac{1}{N^2}\sum_{t,t'} \Cov\!\parens{f(x^{(t)}), f(x^{(t')})}
= \frac{\sigma^2}{N}\parens{1 + 2\sum_{k \ge 1}\parens{1 - \tfrac{k}{N}}\rho(k)}
\;\approx\; \frac{\sigma^2}{N}\,\tau,
$$

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 \abs{\lambda_2}^k$ is geometric, $\tau \approx (1 + \abs{\lambda_2})/(1 - \abs{\lambda_2})$,
which blows up as $\abs{\lambda_2} \to 1$. Thus $N$ MCMC steps are worth only
$N/\tau$ iid samples. A chain **mixes** quickly when $\abs{\lambda_2}$ is small
($\tau$ small) and slowly when $\abs{\lambda_2} \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.

$$
% caption: Slow vs. fast mixing on a bimodal target: left, a trapped chain covers only one mode; right, a well-mixed chain crosses the barrier and visits both.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % ---- LEFT panel: stuck ----
  \begin{scope}
    \draw[->, thick] (-0.2,0) -- (4.6,0) node[right, font=\footnotesize] {$x$};
    \draw[->, thick] (0,-0.2) -- (0,2.6);
    \node[font=\footnotesize, anchor=south] at (2.2,2.5) {slow: stuck in one mode};
    % bimodal density
    \draw[acc, very thick] plot[domain=0.2:4.2, samples=90]
      (\x, {1.7*exp(-(\x-1.2)*(\x-1.2)/0.25) + 1.7*exp(-(\x-3.2)*(\x-3.2)/0.25)});
    % barrier label
    \node[red, font=\scriptsize, anchor=north] at (2.2,-0.28) {barrier};
    \draw[red, thick, <->] (2.0,0.12) -- (2.4,0.12);
    % samples all in left mode
    \foreach \sx in {1.0,1.15,1.3,1.05,1.25,1.18,1.35,1.1,1.22,1.28}
      \fill[green] (\sx,0.07) circle (1.6pt);
  \end{scope}
  % ---- RIGHT panel: mixing ----
  \begin{scope}[xshift=6.0cm]
    \draw[->, thick] (-0.2,0) -- (4.6,0) node[right, font=\footnotesize] {$x$};
    \draw[->, thick] (0,-0.2) -- (0,2.6);
    \node[font=\footnotesize, anchor=south] at (2.2,2.5) {fast: visits both modes};
    \draw[acc, very thick] plot[domain=0.2:4.2, samples=90]
      (\x, {1.7*exp(-(\x-1.2)*(\x-1.2)/0.25) + 1.7*exp(-(\x-3.2)*(\x-3.2)/0.25)});
    % samples in both modes
    \foreach \sx in {1.0,1.2,1.3,1.15,1.28}
      \fill[green] (\sx,0.07) circle (1.6pt);
    \foreach \sx in {3.0,3.2,3.3,3.15,3.28}
      \fill[green] (\sx,0.07) circle (1.6pt);
    \draw[green, thick, ->] (1.5,0.45) .. controls (2.2,1.0) .. (2.9,0.45);
    \node[green, font=\scriptsize, anchor=south] at (2.2,0.85) {crosses};
  \end{scope}
\end{tikzpicture}
$$

| Symptom | Cause | Remedy |
| --- | --- | --- |
| Long burn-in | $\abs{\lambda_2}$ near $1$; poor start | discard early samples; better init |
| High autocorrelation | small proposal steps; strong coupling | thinning; block / reparameterized moves |
| Trapped in one mode | energy barrier between modes | tempering; parallel / annealed chains |
| Slow Gibbs | strongly correlated coordinates | block Gibbs; rotate coordinates |

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.

$$
\nabla_\theta \log p(x)
= \nabla_\theta \log \tilde p(x) - \nabla_\theta \log Z.
$$

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

$$
\nabla_\theta \log Z
= \frac{1}{Z}\nabla_\theta Z
= \frac{1}{Z}\sum_x \nabla_\theta \tilde p(x)
= \sum_x \frac{\tilde p(x)}{Z}\,\nabla_\theta \log \tilde p(x)
= \mathbb{E}_{x \sim p}\!\brackets{\nabla_\theta \log \tilde p(x)}.
$$

> **Theorem (Positive and negative phases).** The log-likelihood gradient of an
> energy-based model is the difference between a data expectation and a model
> expectation,
> $$
> \nabla_\theta \mathbb{E}_{x \sim p_{\text{data}}}[\log p(x)]
> = \underbrace{\mathbb{E}_{x \sim p_{\text{data}}}\!\brackets{\nabla_\theta \log \tilde p(x)}}_{\text{positive phase}}
> - \underbrace{\mathbb{E}_{x \sim p_{\text{model}}}\!\brackets{\nabla_\theta \log \tilde p(x)}}_{\text{negative phase}}.
> $$

> **Proof.** Take $\nabla_\theta \log p(x) = \nabla_\theta \log \tilde p(x) -
> \nabla_\theta \log Z$ and substitute the identity above for the second term, then
> average the first term over $p_{\text{data}}$. The positive phase is an average
> over _data_ — cheap. The negative phase is an average over the _model_ $p$, which
> has no direct sampler — precisely the MCMC problem. $\qed$

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.[^gf-partition]

$$
% caption: The two phases of energy-based learning: the positive phase lowers energy at data, the negative phase raises it at MCMC samples from the model.
\begin{tikzpicture}[>=stealth, font=\small, scale=1.0]
  \definecolor{acc}{HTML}{2348F2}
  \definecolor{green}{HTML}{1F9D4D}
  \definecolor{red}{HTML}{C0392B}
  % axes
  \draw[->, thick] (-0.2,0) -- (6.4,0) node[right, font=\footnotesize] {$x$};
  \draw[->, thick] (0,-0.2) -- (0,3.0) node[above, font=\footnotesize] {energy $E(x)$};
  % energy curve (double well)
  \draw[black, very thick] plot[domain=0.3:6.0, samples=110]
    (\x, {1.4 + 1.1*sin(1.5*\x r) + 0.18*\x});
  % data point: push energy DOWN
  \fill[green] (1.7,1.05) circle (2.6pt);
  \node[green, font=\footnotesize, anchor=south] at (1.7,1.2) {data};
  \draw[->, green, very thick] (1.7,0.95) -- (1.7,0.25);
  \node[green, font=\scriptsize, anchor=west] at (1.8,0.55) {push down};
  % model sample: push energy UP
  \fill[red] (4.5,1.7) circle (2.6pt);
  \node[red, font=\footnotesize, anchor=south] at (3.15,2.15) {model sample};
  \draw[->, red, very thick] (4.5,1.85) -- (4.5,2.65);
  \node[red, font=\scriptsize, anchor=west] at (4.6,2.3) {push up};
  % MCMC label
  \node[font=\scriptsize, black, anchor=west] at (3.2,0.35) {(sampled by MCMC)};
\end{tikzpicture}
$$

### 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 \tfrac1n\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$.

> **Definition (Annealed importance sampling).** A method estimating the ratio
> $Z_T/Z_0$ as a telescoping product of intermediate ratios. A particle is sampled
> from $p_0$ and carried through the bridge by an MCMC transition $T_j$ at each rung;
> the accumulated weight
> $w = \prod_{j=1}^{T} \tilde p_j(x_{j-1}) / \tilde p_{j-1}(x_{j-1})$ is an unbiased
> estimator of $Z_T/Z_0$, so $Z_T = Z_0\,\mathbb{E}[w]$.

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.[^gf-ais]

## 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' = x + \tfrac{\epsilon^2}{2}\nabla_x \log \tilde p(x) + \epsilon\,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 = \tfrac1n\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\,\abs{f}$, 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')g(x\mid x')/(\tilde p(x)g(x'\mid x)))$ — $Z$ cancels;
  **Gibbs** resamples each coordinate from its full conditional, leaving $p$ invariant.
- **Mixing** is governed by $\abs{\lambda_2}$: 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](/deep-learning/generative-models/energy-based-and-boltzmann-machines).

[^gf-mc]: **Goodfellow**, _Deep Learning_, §17.1 — Sampling and Monte Carlo Methods: the sample-mean estimator, its unbiasedness, and the dimension-free $\sigma/\sqrt n$ error rate.
[^gf-importance]: **Goodfellow**, _Deep Learning_, §17.2 — Importance Sampling: the reweighting identity, the optimal proposal $q^\star \propto p\abs{f}$, and the variance blow-up from light-tailed proposals.
[^gf-markov]: **Goodfellow**, _Deep Learning_, §17.3 — Markov Chain Monte Carlo Methods: transition kernels, stationary distributions, and convergence governed by the second eigenvalue.
[^gf-gibbs]: **Goodfellow**, _Deep Learning_, §17.4 — Gibbs Sampling: coordinate-wise resampling from full conditionals, the natural sampler for the sigmoid units of an RBM.
[^gf-partition]: **Goodfellow**, _Deep Learning_, §18.1 — The Log-Likelihood Gradient: the positive/negative phase split and the model-expectation term that demands MCMC.
[^gf-ais]: **Goodfellow**, _Deep Learning_, §18.7 — Annealed Importance Sampling: bridging a tractable base to the target through intermediate distributions to estimate $Z$.