Markov Chain Monte Carlo: Metropolis-Hastings and Gibbs Sampling

Intuition Beginner

Bayesian inference hands you a posterior distribution — a description of how plausible each possible value of an unknown quantity is, given the data. The trouble is that you can almost never write that distribution down cleanly. Bayes' rule gives you the shape of the posterior up to a missing scaling number, and computing that number means doing an integral over every possible parameter value at once. In one or two dimensions you might manage it; in the hundreds of dimensions that real models have, the integral is hopeless. So you are stuck holding a curve whose height you know only up to an unknown overall stretch.

Markov chain Monte Carlo is the escape. Instead of computing the distribution, you take a long random walk through the space of parameter values, arranged so that the walker lingers in regions where the posterior is high and hurries through regions where it is low. After enough steps, the places the walker has visited are spread out in exactly the same proportions as the posterior itself. You never needed the missing scaling number, because the walk only ever compares two heights — it asks "is this new spot more or less plausible than where I am now?" — and the unknown stretch cancels in the comparison.

Picture a politician touring islands of different populations, wanting to spend time on each island in proportion to how many people live there, but without ever being told the actual head counts. Each day they consider hopping to a neighbouring island. If it is more populous, they go. If it is less populous, they flip a biased coin and go only sometimes, with a chance set by the ratio of the two populations. Follow this rule for years and the fraction of nights spent on each island matches its share of the total population, even though no total was ever counted. That ratio-only rule is the heart of every method here.

The one-sentence takeaway: MCMC turns an impossible integral into a long guided random walk, because the long-run share of time the walk spends in each region equals the posterior probability of that region, and only ratios of probabilities are ever needed.

Visual Beginner

The picture below shows a walker exploring a two-humped target and the histogram it builds up over time.

The landscape is the target distribution, drawn at an unknown overall height. The walker spends most of its visits under the tall hump, fewer under the short hump, and rarely crosses the valley — so the histogram of where it has been (bottom) reproduces the landscape's shape (top). The bracket marks the early steps, called burn-in, that are thrown away because the walk had not yet reached the high-probability region.

Worked example Beginner

A coin has unknown success probability $p$. You see $7$ heads in $10$ flips. With a flat prior, the posterior for $p$ is proportional to $p^7(1-p{)}^3$ on the interval from $0$ to $1$. We sample from it with a small Metropolis walk using a fixed random sequence, so you can follow every step.

Step 1. Start at $p=0.5$. The target height there is $0.5^7\times 0.5^3=0.5^{10}=0.000977$.

Step 2. Propose a nearby value by adding a small jump; suppose the proposal is $p^{\prime }=0.7$. The target height there is $0.7^7\times 0.3^3=0.0823\times 0.027=0.002223$.

Step 3. Form the ratio of new height to current height: $0.002223/0.000977=2.28$. The proposal is more plausible than where we are, so the ratio is above $1$. The rule says: when the ratio is at least $1$, always accept. Move to $p=0.7$.

Step 4. From $p=0.7$, propose $p^{\prime }=0.55$. Its height is $0.55^7\times 0.45^3=0.01522\times 0.0911=0.001387$. The ratio is $0.001387/0.002223=0.624$, which is below $1$.

Step 5. When the ratio is below $1$, accept only with probability equal to the ratio. Draw a uniform random number; say it is $0.41$. Since $0.41<0.624$, accept and move to $p=0.55$. Had the random number been $0.80$, we would have stayed at $p=0.7$ and recorded $0.7$ again.

What this tells us: the walk drifted from $0.5$ toward the region around $0.7$, which is where $p^7(1-p{)}^3$ peaks (its maximum is at $p=0.7$). Recording many such steps and making a histogram reproduces the posterior, and at no point did we compute the constant that makes $p^7(1-p{)}^3$ integrate to $1$ — only ratios were used.

Check your understanding Beginner

Formal definition Intermediate+

The target is a probability distribution $\pi$ on a state space $E$ (for Bayesian work, $\pi$ is the posterior $\pi (\theta )\propto L(\theta )p(\theta )$ on a parameter space, known only up to the marginal likelihood $Z=\int L(\theta )p(\theta )d\theta$). Write $\pi (x)$ for the (possibly unnormalised) target density against a reference measure. The aim is to compute expectations $\pi (f)=\int fd\pi$ for test functions $f$, using only evaluations of $\pi$ up to its unknown constant.

Definition (MCMC sampler). A Markov chain Monte Carlo sampler for $\pi$ is a Markov transition kernel $P(x,\cdot )$ on $E$, in the sense of 37.05.01, whose stationary distribution is $\pi$ — that is, $\pi P=\pi$, meaning $\int \pi (dx)P(x,A)=\pi (A)$ for every measurable $A$. The chain $(X_n)$ run from any start is used to form the ergodic average ${\hat{\pi }}_n(f)=\frac{1}{n}{\sum }_{k=0}^{n-1}f(X_k)$, the MCMC estimator of $\pi (f)$.

Definition (reversibility / detailed balance). A kernel $P$ is reversible with respect to $\pi$, or satisfies detailed balance, if

$$\pi (dx)P(x,dy)=\pi (dy)P(y,dx)$$

as measures on $E\times E$. Detailed balance implies $\pi$-invariance: integrating over the first coordinate gives $\int \pi (dx)P(x,dy)=\pi (dy)\int P(y,dx)=\pi (dy)$. It is the standard design tool because it is a pointwise symmetry one can enforce locally, established for Markov chains in 37.05.07.

Definition (Metropolis-Hastings kernel). Fix a proposal kernel $q(x,\cdot )$ with density $q(x,y)$. Given $X_n=x$, draw a proposal $Y\sim q(x,\cdot )$ and accept it with probability

$$\alpha (x,y)=\min \left(1,\frac{\pi (y)q(y,x)}{\pi (x)q(x,y)}\right),$$

setting $X_{n+1}=Y$ on acceptance and $X_{n+1}=x$ on rejection. The resulting kernel is

$$P(x,dy)=q(x,y)\alpha (x,y)dy+r(x){\delta }_x(dy),r(x)=1-\int q(x,y)\alpha (x,y)dy,$$

a mixture of the accepted-move density and a point mass at $x$ of total rejection mass $r(x)$. Only the ratio $\pi (y)/\pi (x)$ enters $\alpha$, so the normalising constant $Z$ never appears.

Two special cases. Random-walk Metropolis takes a symmetric proposal $q(x,y)=q(y,x)$ (for example $Y=x+\epsilon$, $\epsilon$ symmetric), so $\alpha =\min (1,\pi (y)/\pi (x))$ — the original Metropolis form. The independence sampler takes $q(x,y)=g(y)$ independent of $x$, giving $\alpha =\min (1,[\pi (y)g(x)]/[\pi (x)g(y)])$, which is efficient only when $g$ approximates $\pi$ well.

Definition (Gibbs sampler). Suppose $x=(x_1,\dots ,x_d)$ and the full conditionals $\pi (x_i\mid x_{-i})$ — the conditional target of component $i$ given all others — are available to sample. The Gibbs sampler updates components in turn, replacing $x_i$ by a draw from $\pi (\cdot \mid x_{-i})$. The deterministic-scan cycle sweeps $i=1,\dots ,d$; the random-scan version picks $i$ uniformly each step. Each component update is a Metropolis-Hastings step with proposal equal to the full conditional, whose acceptance probability is identically $1$.

Counterexamples to common slips Intermediate+

• Detailed balance is sufficient but not necessary for invariance. A directed cycle on three states with equal stationary weights is $\pi$-invariant yet not reversible. MCMC uses detailed balance because it is constructive, not because invariance demands it.
• A reversible chain need not be ergodic. If the proposal cannot reach part of the support — a random walk with steps too small to cross a deep valley in finite practical time, or zero proposal density on a region — the chain is $\pi$-invariant but explores only one mode. Invariance is about the target; convergence additionally needs $\pi$-irreducibility and aperiodicity.
• The proposal ratio $q(y,x)/q(x,y)$ is not optional for asymmetric proposals. Dropping it (using the Metropolis form $\min (1,\pi (y)/\pi (x))$ with an asymmetric $q$) produces a chain whose stationary law is not $\pi$. The Hastings correction is exactly what restores detailed balance.
• Gibbs updates change one block and condition on the current value of the rest. Conditioning on stale values from before the sweep began, rather than the most recently updated components, breaks invariance of the cycle scan.

Key theorem with proof Intermediate+

Theorem (Metropolis-Hastings invariance). Let $\pi$ be a target with density $\pi (x)$ and $q(x,\cdot )$ a proposal kernel with $q(x,y)>0⟺q(y,x)>0$. The Metropolis-Hastings kernel $P$ with acceptance probability $\alpha (x,y)=\min (1,[\pi (y)q(y,x)]/[\pi (x)q(x,y)])$ is reversible with respect to $\pi$, hence $\pi P=\pi$. The constant $Z$ in $\pi =\tilde{\pi }/Z$ does not affect $P$.

Proof. It suffices to verify detailed balance for the off-diagonal (move) part of the kernel, since the rejection part $r(x){\delta }_x(dy)$ is supported on the diagonal $x=y$, where detailed balance holds automatically (both sides equal $\pi (dx)r(x)$ when $y=x$). For $x\ne y$ the move density of $P$ is $\pi (x)q(x,y)\alpha (x,y)$ against $dxdy$, and detailed balance requires

$$\pi (x)q(x,y)\alpha (x,y)=\pi (y)q(y,x)\alpha (y,x).\text{\hspace{1em}(†)}$$

Fix $x\ne y$ with $q(x,y)>0$, and set $\rho =\frac{\pi (y)q(y,x)}{\pi (x)q(x,y)}$, so $\alpha (x,y)=\min (1,\rho )$. By symmetry of the ratio, $\alpha (y,x)=\min (1,1/\rho )$. Consider the two cases.

If $\rho \le 1$: the left side of $(†)$ is $\pi (x)q(x,y)\cdot \rho =\pi (x)q(x,y)\cdot \frac{\pi (y)q(y,x)}{\pi (x)q(x,y)}=\pi (y)q(y,x)$. The right side is $\pi (y)q(y,x)\cdot \min (1,1/\rho )=\pi (y)q(y,x)\cdot 1$ since $1/\rho \ge 1$. Both sides equal $\pi (y)q(y,x)$.

If $\rho >1$: the left side is $\pi (x)q(x,y)\cdot 1=\pi (x)q(x,y)$. The right side is $\pi (y)q(y,x)\cdot \min (1,1/\rho )=\pi (y)q(y,x)\cdot \frac{1}{\rho }=\pi (y)q(y,x)\cdot \frac{\pi (x)q(x,y)}{\pi (y)q(y,x)}=\pi (x)q(x,y)$. Both sides equal $\pi (x)q(x,y)$.

So $(†)$ holds for all $x\ne y$, and the kernel is reversible with respect to $\pi$. Reversibility gives invariance, $\pi P=\pi$, as in the definition. Finally, replacing $\pi$ by $c\pi$ for any constant $c>0$ leaves $\rho$ — a ratio of $\pi$-values — unchanged, so $P$ depends on $\pi$ only through its shape; the unknown normaliser $Z$ drops out entirely. $\square$

Bridge. This invariance proof builds toward the ergodic theorem that makes ergodic averages consistent, and it appears again in every sampler below, because the Gibbs sampler is the special case in which each proposal is a full conditional and $(†)$ holds with both acceptance probabilities equal to $1$. The foundational reason the algorithm needs only ratios is structural: detailed balance is a statement about the symmetric measure $\pi (dx)P(x,dy)$, and a common rescaling of $\pi$ cancels on the two sides of $(†)$. This is exactly the cancellation foreshadowed in the Beginner island picture, now made into an identity, and it generalises the reversible-chain construction of 37.05.07 from a finite state space to a general one. Putting these together, the design recipe for any MCMC sampler is to engineer detailed balance for the desired target and then invoke ergodicity; the central insight is that local reversibility, checked one pair $(x,y)$ at a time, certifies a global invariant law. The bridge to computation is that this invariant law is the posterior of 26.07.01, so a chain you can simulate replaces an integral you cannot evaluate.

Exercises Intermediate+

Advanced results Master

The invariance theorem gives a chain with the right stationary law; turning ergodic averages into consistent, asymptotically normal estimators requires the general-state-space ergodic theory, and tuning the sampler requires understanding its rate.

Theorem 1 (consistency of ergodic averages). Let $P$ be a $\pi$-invariant transition kernel that is $\pi$-irreducible and Harris recurrent. For every $\pi$-integrable $f$ and every starting point,

$${\hat{\pi }}_n(f)=\frac{1}{n}\sum _{k=0}^{n-1}f(X_k)\stackrel{a.s.}{\to }\pi (f)=\int fd\pi .$$

This is the general-state-space analogue of the Markov-chain ergodic theorem of 37.05.07; Harris recurrence replaces positive recurrence as the hypothesis ensuring that the time-average is independent of the start. For Metropolis-Hastings chains with a proposal positive on the support of $\pi$, $\pi$-irreducibility and Harris recurrence are met under mild conditions (Tierney 1994 ^{[Tierney 1994]}), so the constructed sampler is automatically a valid estimator of posterior expectations.

Theorem 2 (MCMC central limit theorem and the spectral gap). Suppose $P$ is $\pi$-reversible with a spectral gap (its $L^2(\pi )$ self-adjoint operator has spectral radius below $1$ on mean-zero functions), and $f\in L^2(\pi )$. Then

$$\sqrt{n}({\hat{\pi }}_n(f)-\pi (f))\stackrel{d}{\to }\mathcal{N}(0,{\sigma }_f^2),{\sigma }_f^2={\mathrm{Var}}_{\pi }(f)+2\sum _{k=1}^{\infty }{\mathrm{Cov}}_{\pi }(f(X_0),f(X_k)).$$

The asymptotic variance is the integrated autocovariance; for a reversible chain it equals $⟨f,(I+P)(I-P{)}^{-1}f{⟩}_{\pi }$ on mean-zero $f$, finite precisely when the gap is positive. The effective sample size $n_{\mathrm{eff}}=n{\mathrm{Var}}_{\pi }(f)/{\sigma }_f^2=n/{\tau }_f$ with integrated autocorrelation time ${\tau }_f=1+2{\sum }_{k\ge 1}\mathrm{Corr}(f(X_0),f(X_k))$ quantifies how many independent draws the correlated chain is worth.

Theorem 3 (geometric ergodicity and convergence rate). If $P$ is $\pi$-reversible and there exist $V\ge 1$, $\lambda <1$, $b<\infty$ with the drift condition $PV\le \lambda V+b{1}_C$ for a small set $C$, then $P$ is geometrically ergodic: $\Vert P^n(x,\cdot )-\pi {\Vert }_{TV}\le M(x){\rho }^n$ for some $\rho <1$. Geometric ergodicity is what licences a finite burn-in and the MCMC central limit theorem for all $f$ with $\pi (f^{2}V)<\infty$. Random-walk Metropolis on a target with exponentially decaying tails is geometrically ergodic; on heavy-tailed targets it can fail to be, and the ergodic average then converges with no Gaussian-rate guarantee.

Theorem 4 (optimal scaling of random-walk Metropolis). For a product target in dimension $d$ with proposal $Y=x+{\sigma }_{d}Z$, $Z$ standard, the efficient scaling is ${\sigma }_d\propto d^{-1/2}$, and as $d\to \infty$ the rescaled first coordinate converges to a Langevin diffusion whose speed is maximised at an acceptance rate of $0.234$ (Roberts-Gelman-Gilks 1997 ^{[Roberts-Gelman-Gilks 1997]}). The practical tuning heuristic — aim for roughly a quarter of proposals accepted — is the finite-dimensional shadow of this diffusion limit.

Synthesis. The foundational reason MCMC works is a clean separation of two problems that detailed balance and ergodicity solve in turn: invariance is engineered locally by the Metropolis-Hastings acceptance rule, and consistency is then inherited globally from the ergodic theorem, exactly as in the finite-chain theory of 37.05.07 now lifted to a continuum by Harris recurrence. This is the central insight of the subject — the same trajectory that explores the posterior also estimates expectations over it, so one simulation does the work of both sampling and integration. The asymptotic variance of Theorem 2 and the geometric rate of Theorem 3 are dual descriptions of a single quantity, the spectral gap: the gap controls how fast $P^n(x,\cdot )$ forgets its start and, identically, how small the integrated autocorrelation time of the estimator is. Putting these together with the optimal-scaling result of Theorem 4 gives the complete operational picture, where tuning the proposal is tuning the gap and the gap is what the $0.234$ acceptance heuristic silently optimises. The bridge from this machinery back to inference is that the invariant law is the Bayesian posterior of 26.07.01, and the ratio-only acceptance rule generalises the reversible-chain construction so that the intractable marginal likelihood $Z$ never has to be computed; what looks like an integration problem is, on inspection, a simulation problem.

Full proof set Master

Proposition 1 (detailed balance implies invariance, general state space). If $\pi (dx)P(x,dy)=\pi (dy)P(y,dx)$ as measures on $E^2$, then $\pi P=\pi$.

Proof. For measurable $A\subseteq E$, integrate the symmetric measure over $\{x\in E,y\in A\}$. The left side gives ${\int }_E\pi (dx)P(x,A)=(\pi P)(A)$. The right side gives ${\int }_{y\in A}\pi (dy){\int }_{x\in E}P(y,dx)={\int }_A\pi (dy)\cdot 1=\pi (A)$, since $P(y,\cdot )$ is a probability measure. Hence $(\pi P)(A)=\pi (A)$ for all $A$. $\square$

Proposition 2 (Metropolis-Hastings reversibility). The kernel of the Formal-definition section satisfies detailed balance for $\pi$.

Proof. On the diagonal the rejection mass $r(x){\delta }_x$ contributes symmetrically. Off-diagonal, the claim is $\pi (x)q(x,y)\alpha (x,y)=\pi (y)q(y,x)\alpha (y,x)$ for $x\ne y$. With $\rho =[\pi (y)q(y,x)]/[\pi (x)q(x,y)]$ and $\alpha (x,y)=\min (1,\rho )$, $\alpha (y,x)=\min (1,1/\rho )$: if $\rho \le 1$ both sides equal $\pi (y)q(y,x)$; if $\rho >1$ both sides equal $\pi (x)q(x,y)$. The two cases are the computation in the Key theorem. Reversibility plus Proposition 1 give invariance. $\square$

Proposition 3 (Gibbs as rejection-free Metropolis-Hastings, with invariance). A full-conditional update of component $i$, $x\mapsto (x_i^{\prime },x_{-i})$ with $x_i^{\prime }\sim \pi (\cdot \mid x_{-i})$, is a Metropolis-Hastings step with acceptance probability $1$, and it leaves $\pi$ invariant.

Proof. Take proposal $q(x,x^{\prime })=\pi (x_i^{\prime }\mid x_{-i})$ for moves changing only coordinate $i$. The acceptance ratio is

$$\frac{\pi (x^{\prime })q(x^{\prime },x)}{\pi (x)q(x,x^{\prime })}=\frac{\pi (x_i^{\prime }\mid x_{-i})\pi (x_{-i})\cdot \pi (x_i\mid x_{-i})}{\pi (x_i\mid x_{-i})\pi (x_{-i})\cdot \pi (x_i^{\prime }\mid x_{-i})}=1,$$

using $\pi (x)=\pi (x_i\mid x_{-i})\pi (x_{-i})$ and that $x,x^{\prime }$ share $x_{-i}$. So $\alpha \equiv 1$. Invariance follows from Proposition 2, or directly: the kernel $P_i(x,d{x}^{\prime })=\pi (x_i^{\prime }\mid x_{-i}){\delta }_{x_{-i}}(d{x}_{-i}^{\prime })$ has $\int \pi (dx)P_i(x,A)=\int \pi (x_{-i})[\int \pi (x_i\mid x_{-i})d{x}_i]\pi (x_i^{\prime }\mid x_{-i})1_{A}d{x}_i^{\prime }d{x}_{-i}=\int \pi (x_i^{\prime }\mid x_{-i})\pi (x_{-i})1_A=\pi (A)$, the inner $x_i$-integral being $1$. $\square$

Proposition 4 (invariance of compositions and mixtures). If $P_1,\dots ,P_d$ are each $\pi$-invariant, then both the composition $P_1\cdots P_d$ and any convex mixture ${\sum }_j{\beta }_j{P}_j$ are $\pi$-invariant.

Proof. Composition: $\pi (P_1\cdots P_d)=(\pi P_1)P_2\cdots P_d=\pi P_2\cdots P_d=\cdots =\pi$. Mixture: $\pi \left({\sum }_j{\beta }_j{P}_j\right)={\sum }_j{\beta }_j(\pi P_j)={\sum }_j{\beta }_j\pi =\pi$, since ${\sum }_j{\beta }_j=1$. The cycle-scan Gibbs sampler is the composition $P_1\cdots P_d$; the random-scan version is the mixture $\frac{1}{d}{\sum }_i{P}_i$. Both are $\pi$-invariant by Proposition 3 plus this proposition. $\square$

Proposition 5 (normalising constant cancels). Replacing the target $\pi$ by $c\pi$ for any $c>0$ leaves the Metropolis-Hastings kernel unchanged.

Proof. The kernel depends on $\pi$ only through $\alpha (x,y)=\min (1,[\pi (y)q(y,x)]/[\pi (x)q(x,y)])$. The substitution $\pi \mapsto c\pi$ multiplies numerator and denominator of the inner ratio by $c$, which cancels, so $\alpha$ and hence $P$ are unchanged. In particular, knowing the posterior only up to its marginal likelihood $Z$ is sufficient to run the sampler. $\square$

Proposition 6 (asymptotic variance for a reversible chain). For a $\pi$-reversible $P$ with spectral gap and mean-zero $f\in L^2(\pi )$, the limiting variance of $\sqrt{n}{\hat{\pi }}_n(f)$ is ${\sigma }_f^2=⟨f,(I+P)(I-P{)}^{-1}f{⟩}_{\pi }$, which is nonnegative and finite.

Proof. By the spectral theorem for the self-adjoint operator $P$ on mean-zero $L^2(\pi )$, write $f=\int \varphi d{E}_{\lambda }f$ over the spectrum $\lambda \in (-1,1)$ (gap excludes the endpoints from the mean-zero part). Stationary lag covariances are $\mathrm{Cov}(f(X_0),f(X_k))=⟨f,P^{k}f⟩=\int {\lambda }^{k}d{\mu }_f(\lambda )$ with ${\mu }_f$ the spectral measure of $f$. Summing the geometric series, ${\sigma }_f^2=\int (1+2{\sum }_{k\ge 1}{\lambda }^k)d{\mu }_f=\int \frac{1+\lambda }{1-\lambda }d{\mu }_f(\lambda )=⟨f,(I+P)(I-P{)}^{-1}f⟩$. The integrand $\frac{1+\lambda }{1-\lambda }\ge 0$ on $(-1,1)$ and is bounded by $\frac{1+{\lambda }_{\max }}{1-{\lambda }_{\max }}<\infty$ under the gap, so ${\sigma }_f^2$ is finite and nonnegative. The finite-state eigen-sum of Exercise 8 is the discrete case. $\square$

Connections Master

• The Markov property and transition matrices of 37.05.01 supply the kernel formalism on which every sampler is built: an MCMC sampler is a transition kernel, and the Chapman-Kolmogorov composition there is what lets the cycle-scan Gibbs sampler be analysed as a product of component kernels. The $n$-step convergence of that unit is the finite-state shadow of the geometric ergodicity used here to license burn-in.

• The ergodic theorem and detailed balance of 37.05.07 are the direct parent: that unit proves, on a countable state space, both that detailed balance gives a stationary law and that time-averages converge to space-averages, and it already exhibits the Metropolis-Hastings acceptance ratio as the reversibility-design tool. The present unit specialises the construction to a target posterior and lifts the convergence statement to a general state space via Harris recurrence, so the algorithm estimates expectations rather than merely admitting a stationary law.

• The Bayesian prior-posterior framework of 26.07.01 is what MCMC makes computational: the target $\pi$ here is precisely the posterior $\pi (\theta )\propto L(\theta )p(\theta )$ defined there, and the cancellation of the normalising constant is exactly what removes the intractable marginal likelihood from posterior summarisation, turning conjugate-only Bayesian analysis into a method usable on any model whose likelihood can be evaluated.

Historical & philosophical context Master

The algorithm originates in statistical physics. Metropolis, the Rosenbluths, and the Tellers introduced the symmetric-proposal accept/reject scheme in 1953 to compute equilibrium properties of interacting particles by sampling the Boltzmann distribution on the MANIAC computer ^{[Metropolis 1953]}, formulating the rule as an acceptance probability depending only on the energy difference — the physical analogue of the target ratio. Hastings generalised the method to arbitrary, possibly asymmetric, proposals in 1970, supplying the correction factor $q(y,x)/q(x,y)$ that makes the construction work for any target and recasting it as a tool for general Monte Carlo computation rather than a physics-specific device ^{[Hastings 1970]}.

The method entered statistics decisively through image analysis: Geman and Geman introduced the Gibbs sampler in 1984 for Bayesian restoration of images modelled by Gibbs random fields, naming it for the Gibbs distributions it sampled and exploiting the tractable full conditionals of Markov random fields ^{[Geman 1984]}. The recognition that the same machinery applied to general Bayesian posteriors, made widely known by Gelfand and Smith in 1990, set off the computational revolution that made hierarchical and high-dimensional Bayesian models practical. The general-state-space theory — $\pi$-irreducibility, Harris recurrence, and the ergodic theorem and central limit theorem that justify the estimators — was systematised by Tierney in 1994 ^{[Tierney 1994]}, building on the Markov-chain stability theory of Meyn and Tweedie.

Bibliography Master

@article{Metropolis1953,
  author  = {Metropolis, Nicholas and Rosenbluth, Arianna W. and Rosenbluth, Marshall N. and Teller, Augusta H. and Teller, Edward},
  title   = {Equation of State Calculations by Fast Computing Machines},
  journal = {Journal of Chemical Physics},
  volume  = {21},
  number  = {6},
  pages   = {1087--1092},
  year    = {1953}
}

@article{Hastings1970,
  author  = {Hastings, W. K.},
  title   = {Monte Carlo Sampling Methods Using Markov Chains and Their Applications},
  journal = {Biometrika},
  volume  = {57},
  number  = {1},
  pages   = {97--109},
  year    = {1970}
}

@article{GemanGeman1984,
  author  = {Geman, Stuart and Geman, Donald},
  title   = {Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images},
  journal = {IEEE Transactions on Pattern Analysis and Machine Intelligence},
  volume  = {6},
  number  = {6},
  pages   = {721--741},
  year    = {1984}
}

@article{GelfandSmith1990,
  author  = {Gelfand, Alan E. and Smith, Adrian F. M.},
  title   = {Sampling-Based Approaches to Calculating Marginal Densities},
  journal = {Journal of the American Statistical Association},
  volume  = {85},
  number  = {410},
  pages   = {398--409},
  year    = {1990}
}

@article{Tierney1994,
  author  = {Tierney, Luke},
  title   = {Markov Chains for Exploring Posterior Distributions},
  journal = {Annals of Statistics},
  volume  = {22},
  number  = {4},
  pages   = {1701--1762},
  year    = {1994}
}

@article{RobertsGelmanGilks1997,
  author  = {Roberts, Gareth O. and Gelman, Andrew and Gilks, Walter R.},
  title   = {Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms},
  journal = {Annals of Applied Probability},
  volume  = {7},
  number  = {1},
  pages   = {110--120},
  year    = {1997}
}

@book{RobertCasella2004,
  author    = {Robert, Christian P. and Casella, George},
  title     = {Monte Carlo Statistical Methods},
  edition   = {2},
  publisher = {Springer},
  year      = {2004}
}

@book{GelmanBDA2013,
  author    = {Gelman, Andrew and Carlin, John B. and Stern, Hal S. and Dunson, David B. and Vehtari, Aki and Rubin, Donald B.},
  title     = {Bayesian Data Analysis},
  edition   = {3},
  publisher = {CRC Press},
  year      = {2013}
}

@book{MeynTweedie2009,
  author    = {Meyn, Sean and Tweedie, Richard L.},
  title     = {Markov Chains and Stochastic Stability},
  edition   = {2},
  publisher = {Cambridge University Press},
  year      = {2009}
}