Sampling Methods

When an integral over a posterior distribution cannot be done exactly, we can draw random samples from the distribution and average over them. This chapter builds the samplers: transformations of uniform numbers, rejection and importance sampling, and then Markov chains (Metropolis-Hastings, Gibbs, slice sampling and hybrid Monte Carlo) whose long-run distribution is the one we want. Every proof that the book sketches or leaves to the reader is written out below.

TrungTin Nguyen · Following Bishop (2006), Pattern Recognition and Machine Learning, Chapter 11 · Course home · Mathematical toolkit · Python code for the figures

Ch. 1 Introduction Ch. 2 Probability distributions Ch. 3 Linear regression Ch. 4 Linear classification Ch. 5 Neural networks Ch. 6 Kernel methods Ch. 7 Sparse kernel machines Ch. 8 Graphical models Ch. 9 Mixture models and EM Ch. 10 Approximate inference Ch. 11 Sampling methods Ch. 12 Continuous latent variables Ch. 13 Sequential data Ch. 14 Combining models

Monte Carlo estimates

pp. 523–526 Most of what we want from a posterior distribution is an expectation: a posterior mean, a predictive probability, a posterior probability of a set. So the basic problem of the chapter is to compute

\[ \E[f]=\int f(\mathbf{z})\,p(\mathbf{z})\,\d\mathbf{z} \tag{11.1} \]

for a function \(f\) and a distribution \(p\) (a sum replaces the integral for discrete \(\mathbf{z}\)), when the integral cannot be done analytically. The idea of Monte Carlo is to replace the distribution by a sample: if \(\mathbf{z}^{(1)},\dots,\mathbf{z}^{(L)}\) are drawn independently from \(p\), use

\[ \widehat f=\frac1L\sum_{l=1}^{L} f(\mathbf{z}^{(l)}) . \tag{11.2} \]

Proof

Write \(f_l=f(\mathbf{z}^{(l)})\). Each \(f_l\) has the distribution of \(f(\mathbf{z})\) with \(\mathbf{z}\sim p\), so \(\E[f_l]=\E[f]\) and \(\var[f_l]=\var[f]\).

Mean. Expectation is linear: \(\E[\widehat f]\)\({}=\frac1L\sum_l\E[f_l]\)\({}=\frac1L\cdot L\,\E[f]=\E[f]\). This needs only that each \(\mathbf{z}^{(l)}\) has distribution \(p\), not independence.

Variance. Let \(g_l=f_l-\E[f]\), so \(\E[g_l]=0\) and \(\widehat f-\E[f]=\frac1L\sum_l g_l\). Then

\[ \var[\widehat f]=\E\Big[\Big(\frac1L\sum_{l}g_l\Big)^2\Big]=\frac1{L^2}\sum_{l=1}^{L}\sum_{m=1}^{L}\E[g_lg_m]. \]

For \(l\ne m\), independence gives \(\E[g_lg_m]=\E[g_l]\,\E[g_m]=0\). Only the \(L\) diagonal terms survive, each equal to \(\E[g_l^2]=\E[(f-\E[f])^2]\), so \(\var[\widehat f]\)\({}=\frac{1}{L^2}\cdot L\,\E[(f-\E[f])^2]\), which is (11.3). (Only uncorrelatedness of the \(f_l\) was used.)

Dimension. \(\var[f]\) is a property of the scalar random variable \(f(\mathbf{z})\); the dimension of \(\mathbf{z}\) enters only through it. So the standard error \(\sqrt{\var[f]/L}\) falls like \(L^{-1/2}\) in any dimension, whereas a grid with \(m\) points per axis needs \(m^D\) evaluations. Ten or twenty independent draws give a rough estimate (standard error about a quarter to a third of the standard deviation of \(f\)), which is the sense of the book's remark on p. 524. The limit theorems are the classical ones for independent, identically distributed \(f_l\) with finite mean (strong law) and finite variance (central limit theorem). ∎

The book adds two warnings p. 524. First, the samplers of Section 11.2 produce dependent draws, and then (11.3) no longer holds. Second, if \(f\) is large where \(p\) is small (Figure 11.1), \(\var[f]\) is large and so is the required \(L\). The next result quantifies the first warning; it is the reason why "effective sample size" appears throughout the chapter.

Proof

Unbiasedness was shown above without independence. With \(g_l=f_l-\E[f]\) and stationarity, \(\E[g_lg_m]=\var[f]\,\rho_{|l-m|}\). In the double sum \(\sum_{l,m}\E[g_lg_m]\) the lag \(k=|l-m|\) occurs \(L\) times for \(k=0\) and \(2(L-k)\) times for \(1\le k\le L-1\) (ordered pairs with \(m=l\pm k\)). Hence

\[ \var[\widehat f]=\frac{\var[f]}{L^2}\Big[L+2\sum_{k=1}^{L-1}(L-k)\rho_k\Big], \]

which is the formula. If \(\sum_k|\rho_k|\lt\infty\), then \(\sum_{k\lt L}(1-k/L)\rho_k\to\sum_{k\ge1}\rho_k\) by dominated convergence (each term converges and is bounded by \(|\rho_k|\)). For \(\rho_k=\phi^k\) with \(|\phi|\lt1\): \(1+2\sum_{k\ge1}\phi^k\)\({}=1+2\phi/(1-\phi)\)\({}=(1+\phi)/(1-\phi)\). ∎

Negative correlations give \(\tau_{\rm int}\lt1\): antithetic schemes such as over-relaxation (Section 11.3) can beat independent sampling for some functions. Thinning (keeping every \(M\)th draw, as suggested on p. 538) makes the kept draws nearly independent, but it never reduces the variance of the estimate from a fixed run; it only saves storage.

Ancestral and logic sampling pp. 524–525, (11.4)

For a directed graphical model the joint distribution factorises as

\[ p(\mathbf{z})=\prod_{i=1}^{M}p(\mathbf{z}_i\given \mathrm{pa}_i). \tag{11.4} \]

Ancestral sampling visits the nodes in an order in which parents come before children and draws each \(\mathbf{z}_i\) from \(p(\mathbf{z}_i\given\mathrm{pa}_i)\) with the parents set to their sampled values. The probability of producing \(\mathbf{z}\) is the product of the factors drawn, that is (11.4) itself, so one pass gives an exact draw from the joint. Logic sampling (Henrion, 1988) deals with observed nodes by sampling them too and throwing the whole draw away if a sampled value disagrees with the observation.

For undirected graphs there is no ordering that makes one pass exact; Gibbs sampling (Section 11.3) is used instead. To sample a marginal \(p(\mathbf{u})\) from a joint sampler for \(p(\mathbf{u},\mathbf{v})\), simply ignore \(\mathbf{v}\) p. 525: the \(\mathbf{u}\)-component of a draw from the joint has the marginal distribution by definition.

11.1 Basic sampling algorithms

p. 526 We assume a supply of independent numbers uniform on \((0,1)\). (In practice these are pseudo-random numbers: a deterministic sequence that passes statistical tests of randomness. Their construction is outside the scope of the book and of these notes.) Everything in Section 11.1 transforms, filters or reweights such numbers.

11.1.1 Standard distributions

pp. 526–528 If \(z\) is uniform on \((0,1)\) and \(y=f(z)\) for a smooth monotone \(f\), the density of \(y\) is

\[ p(y)=p(z)\left|\frac{\d z}{\d y}\right| , \tag{11.5} \]

with \(p(z)=1\) on \((0,1)\). We want to choose \(f\) so that \(y\) has a prescribed density \(p(y)\). Integrating (11.5) suggests

\[ z=h(y)\equiv\int_{-\infty}^{y}p(\widehat y)\,\d\widehat y , \tag{11.6} \]

the distribution function of \(y\), and \(y=h^{-1}(z)\) (Figure 11.2).

Erratum. Some printings omit the absolute value in (11.5). It is needed when \(f\) is decreasing.

Proof

The key equivalence. For \(z\in(0,1)\) and \(a\in\R\): \(h^{-}(z)\le a\iff z\le h(a)\). Indeed, if \(z\le h(a)\), then \(a\) belongs to the set \(\{y:h(y)\ge z\}\), so its infimum is at most \(a\). Conversely, distribution functions are right-continuous and nondecreasing, so the set \(\{y:h(y)\ge z\}\) is a closed half-line \([h^-(z),\infty)\); if \(h^{-}(z)\le a\), then \(a\) is in that half-line and \(h(a)\ge z\).

Conclusion. \(\Pr(y\le a)\)\({}=\Pr(z\le h(a))\)\({}=h(a)\), because \(z\) is uniform and \(h(a)\in[0,1]\). So \(y\) has distribution function \(h\). If \(h\) is differentiable with \(h'=p\), the density of \(y\) is \(p\). This is also (11.5) read backwards: with \(z=h(y)\), \(\left|\d z/\d y\right|\)\({}=h'(y)=p(y)\) and \(p(z)=1\). The generalised inverse matters only when \(h\) has flat pieces (gaps in the support) or jumps (discrete distributions); for a discrete distribution on \(y_1\lt y_2\lt\dots\) the rule says: return the first \(y_j\) whose cumulative probability reaches \(z\). ∎

Several variables p. 527. If \(\mathbf{y}=\mathbf{f}(\mathbf{z})\) is a smooth bijection between open sets of \(\R^M\), then

\[ p(y_1,\dots,y_M)=p(z_1,\dots,z_M)\left|\frac{\partial(z_1,\dots,z_M)}{\partial(y_1,\dots,y_M)}\right| , \tag{11.9} \]

where the last factor is the absolute value of the determinant of the matrix of partial derivatives \(\partial z_i/\partial y_j\) (the Jacobian of the inverse map).

The Box-Muller method pp. 527–528, (11.10)–(11.12)

Draw \(z_1,z_2\) uniform on \((-1,1)\) (from \(z\mapsto2z-1\)) and keep the pair only if \(z_1^2+z_2^2\le1\). The kept pairs are uniform on the unit disc, with density \(1/\pi\) (Figure 11.3). With \(r^2=z_1^2+z_2^2\), set

\[ y_1=z_1\left(\frac{-2\ln r^2}{r^2}\right)^{1/2}, \tag{11.10} \] \[ y_2=z_2\left(\frac{-2\ln r^2}{r^2}\right)^{1/2}. \tag{11.11} \]

Then

\[ p(y_1,y_2)=p(z_1,z_2)\left|\frac{\partial(z_1,z_2)}{\partial(y_1,y_2)}\right|=\left[\frac{1}{\sqrt{2\pi}}\exp(-y_1^2/2)\right]\left[\frac{1}{\sqrt{2\pi}}\exp(-y_2^2/2)\right], \tag{11.12} \]

so \(y_1\) and \(y_2\) are independent standard Gaussians.

Erratum. Some printings have \(\ln z_1\) and \(\ln z_2\) in the numerators of (11.10) and (11.11); both should be \(\ln r^2\), as above.

Proof

(a) \((z_1,z_2)\) is uniform on the square \((-1,1)^2\) of area 4. Conditioning a uniform distribution on a subset gives the uniform distribution on that subset, so the kept pairs are uniform on the disc of area \(\pi\), with density \(1/\pi\); the probability of keeping a pair is the ratio of areas, \(\pi/4\). (This is rejection sampling with a uniform proposal, Section 11.1.2.)

(b) Write \(s=r^2=z_1^2+z_2^2\) and \(g(s)=\sqrt{-2\ln s/s}\), so \(y_i=z_i\,g(s)\). The map keeps the direction of \((z_1,z_2)\) and sends the radius \(r\) to \(\rho=r\,g(r^2)=\sqrt{-2\ln r^2}\), which decreases strictly from \(+\infty\) (as \(r\to0\)) to \(0\) (as \(r\to1\)). Hence each ray of the punctured disc is mapped one-to-one onto the corresponding open ray of the plane, and the map is a bijection onto \(\R^2\setminus\{\mathbf 0\}\). For the Jacobian, \(\partial s/\partial z_j=2z_j\), so

\[ \frac{\partial y_i}{\partial z_j}=\delta_{ij}\,g(s)+2z_iz_j\,g'(s),\qquad\text{i.e.}\qquad \mathbf{J}=g\,\mathbf{I}+2g'\,\mathbf{z}\mathbf{z}^\T . \]

The matrix \(\mathbf{z}\mathbf{z}^\T\) has eigenvalue \(s\) (eigenvector \(\mathbf{z}\)) and \(0\) (the perpendicular direction), so \(\mathbf{J}\) has eigenvalues \(g+2sg'\) and \(g\), and \(\det\mathbf{J}=g\,(g+2sg')\). Now \(\ln g=\tfrac12\ln(-2\ln s)-\tfrac12\ln s\), so

\[ \frac{g'}{g}=\frac12\cdot\frac{-2/s}{-2\ln s}-\frac{1}{2s}=\frac{1}{2s\ln s}-\frac1{2s},\qquad g+2sg'=g\Big(1+\frac{1}{\ln s}-1\Big)=\frac{g}{\ln s}. \]

Therefore \(\det\mathbf{J}=g^2/\ln s\)\({}=\frac{-2\ln s}{s}\cdot\frac{1}{\ln s}\)\({}=-\frac{2}{s}\). (A symbolic computation gives the same result.)

(c) By (11.9) with the inverse Jacobian \(|\partial(z_1,z_2)/\partial(y_1,y_2)|\)\({}=1/|\det\mathbf{J}|=s/2\): \(p(y_1,y_2)=\frac1\pi\cdot\frac s2\). It remains to write \(s\) in terms of \(\mathbf{y}\): \(y_1^2+y_2^2=s\,g(s)^2=-2\ln s\), so \(s=\exp\{-(y_1^2+y_2^2)/2\}\) and

\[ p(y_1,y_2)=\frac{1}{2\pi}\exp\Big\{-\frac{y_1^2+y_2^2}{2}\Big\}, \]

which factorises as (11.12). The single point \(\mathbf{y}=\mathbf 0\) and the null event \(s=0\) do not matter.

(d) In polar coordinates the uniform point on the disc has \(s=r^2\) uniform on \((0,1)\) (because \(\Pr(r^2\le t)=\pi t/\pi=t\)) and angle \(\theta\) uniform on \((0,2\pi)\), independent of \(s\). The output has radius \(\sqrt{-2\ln s}\) and the same angle. Replacing \((s,\theta)\) by \((u_1,2\pi u_2)\) gives the trigonometric form. As a by-product, \(y_1^2+y_2^2=-2\ln s\) is exponential with mean 2, the \(\chi^2_2\) distribution, as it must be. The polar form avoids evaluating sine and cosine at the price of rejecting \(1-\pi/4\approx21\%\) of the pairs (Marsaglia and Bray, 1964). ∎

Correlated Gaussians p. 528. If \(y\sim\N(0,1)\), then \(\sigma y+\mu\sim\N(\mu,\sigma^2)\). For a vector with mean \(\bm\mu\) and covariance \(\bm\Sigma\), factorise \(\bm\Sigma=\mathbf{L}\mathbf{L}^\T\) with \(\mathbf{L}\) lower triangular (the Cholesky decomposition) and use \(\mathbf{y}=\bm\mu+\mathbf{L}\mathbf{z}\) with \(\mathbf{z}\) a vector of independent standard Gaussians.

Proof

Moments. \(\E[\mathbf{y}]\)\({}=\bm\mu+\mathbf{A}\E[\mathbf{z}]\)\({}=\bm\mu\), and

\[ \cov[\mathbf{y}]=\E\big[(\mathbf{A}\mathbf{z})(\mathbf{A}\mathbf{z})^\T\big]=\mathbf{A}\,\E[\mathbf{z}\mathbf{z}^\T]\,\mathbf{A}^\T=\mathbf{A}\mathbf{I}\mathbf{A}^\T=\bm\Sigma . \]

Gaussianity. \(\mathbf{A}\) is invertible because \(\det(\mathbf{A})^2=\det\bm\Sigma\gt0\). The inverse map is \(\mathbf{z}=\mathbf{A}^{-1}(\mathbf{y}-\bm\mu)\) with constant Jacobian determinant \(\det\mathbf{A}^{-1}\), so by (11.9)

\[ p(\mathbf{y})=\frac{1}{(2\pi)^{D/2}}\exp\Big\{-\tfrac12(\mathbf{y}-\bm\mu)^\T\mathbf{A}^{-\T}\mathbf{A}^{-1}(\mathbf{y}-\bm\mu)\Big\}\frac{1}{|\det\mathbf{A}|}=\N(\mathbf{y}\given\bm\mu,\bm\Sigma), \]

since \(\mathbf{A}^{-\T}\mathbf{A}^{-1}\)\({}=(\mathbf{A}\mathbf{A}^\T)^{-1}\)\({}=\bm\Sigma^{-1}\) and \(|\det\mathbf{A}|=|\bm\Sigma|^{1/2}\). (See also the toolkit.) Cholesky is the usual choice because it costs \(D^3/3\) operations, is numerically stable for positive definite matrices, and \(\mathbf{L}\mathbf{z}\) costs only \(D^2/2\) multiplications per draw. If the precision \(\bm\Lambda\)\({}=\bm\Sigma^{-1}\)\({}=\mathbf{R}\mathbf{R}^\T\) is given instead, solve \(\mathbf{R}^\T\mathbf{x}=\mathbf{z}\) by back substitution and set \(\mathbf{y}=\bm\mu+\mathbf{x}\): then \(\cov[\mathbf{x}]\)\({}=\mathbf{R}^{-\T}\mathbf{R}^{-1}\)\({}=\bm\Lambda^{-1}\). ∎

The transformation method needs a distribution function that can be computed and inverted, which is rare. The next two techniques need much less: the ability to evaluate \(\widetilde p\) and to sample a simpler \(q\) p. 528.

11.1.2 Rejection sampling

pp. 528–530 Suppose we can evaluate

\[ p(z)=\frac{1}{Z_p}\,\widetilde p(z) \tag{11.13} \]

up to the unknown constant \(Z_p\), and can sample a proposal density \(q(z)\). Choose a constant \(k\) with \(kq(z)\ge\widetilde p(z)\) for all \(z\); the function \(kq\) is the comparison function or envelope (Figure 11.4).

Proof

The area of \(A_g\) is \(\int g(z)\,\d z\), so the uniform density on \(A_g\) is \(\mathbb{1}\{0\le u\le g(z)\}/\int g\). Integrating out \(u\) gives the marginal \(g(z)/\int g\); dividing the joint by this marginal gives the conditional \(\mathbb 1\{0\le u\le g(z)\}/g(z)\), the uniform density on \([0,g(z)]\). Conversely, multiplying that marginal and that conditional gives back the uniform joint density. ∎

Proof

(a) Given \(z_0\), \(u_0\) is uniform on \([0,kq(z_0)]\) and \(0\le\widetilde p(z_0)\le kq(z_0)\), so \(\Pr(u_0\le\widetilde p(z_0)\given z_0)\)\({}=\widetilde p(z_0)/kq(z_0)\). (Points with \(q(z_0)=0\) are never proposed; the envelope condition forces \(\widetilde p=0\) there.)

(b) Average (a) over \(z_0\sim q\): the factors \(q\) cancel, which is (11.14).

(c) By the product rule, the joint density of "propose \(z\) and accept it" is \(q(z)\,\widetilde p(z)/kq(z)\)\({}=\widetilde p(z)/k\). By Bayes' rule the density of \(z\) given acceptance is this divided by (11.14):

\[ p(z\given\text{accept})=\frac{\widetilde p(z)/k}{Z_p/k}=\frac{\widetilde p(z)}{Z_p}=p(z). \]

Geometrically: \((z_0,u_0)\) is uniform under the graph of \(kq\) (Result 11.6 with \(g=kq\)); conditioning on the subset \(A_{\widetilde p}\) leaves it uniform there, and then Result 11.6 with \(g=\widetilde p\) says that \(z_0\) has density \(p\).

(d) Attempts are independent and each succeeds with probability \(Z_p/k\), so the number of attempts is geometric with that success probability and mean \(k/Z_p\). Because each attempt uses fresh random numbers, the accepted value of one round is independent of the others. ∎

So the waste is the fraction of the area under \(kq\) that is not under \(\widetilde p\), and \(k\) should be as small as the condition \(kq\ge\widetilde p\) allows: \(k^\star=\sup_z\widetilde p(z)/q(z)\). If \(p\) has heavier tails than \(q\), the ratio is unbounded and no \(k\) exists.

Erratum. In some printings the label of the lower curve of Figure 11.4 (and of Figure 11.13) is \(p\); it should be \(\widetilde p\).

The gamma distribution with a Cauchy envelope pp. 529–530, (11.15)–(11.16)

The gamma density

\[ \Gam(z\given a,b)=\frac{b^a z^{a-1}\exp(-bz)}{\Gamma(a)} \tag{11.15} \]

is bell-shaped for \(a\gt1\), with mode \((a-1)/b\) (toolkit). A Cauchy density is also bell-shaped, easy to sample by the transformation method, and has heavier tails than the gamma, so a scaled Cauchy can serve as an envelope:

\[ q(z)=\frac{k}{1+(z-c)^2/b^2}, \tag{11.16} \]

which is generated by \(z=b\tan y+c\) with \(y\) uniform on \((-\pi/2,\pi/2)\) (Exercise 11.7). The book states that the rejection rate is smallest for \(c=a-1\), \(b^2=2a-1\), with \(k\) as small as possible. Here the letter \(b\) in (11.16) is the Cauchy scale; for these formulas to make sense the gamma must have rate 1, and we take \(b=1\) in (11.15) for this example (a general rate only rescales \(z\)). The book's Figure 11.5 uses \(a=10\).

11.1.3 Adaptive rejection sampling

pp. 530–532 Finding a tight envelope by hand is often hard. Gilks and Wild (1992) build it automatically when \(p\) is log-concave, that is, when \(\ell(z)=\ln\widetilde p(z)\) is concave (for differentiable \(\ell\): \(\ell'\) is nonincreasing). Many densities met in Gibbs sampling are log-concave: Gaussian, gamma with shape \(\ge1\), beta with both parameters \(\ge1\), logistic regression posteriors, and products of these.

Proof

(a) A differentiable concave function lies below each of its tangents: for concave \(\ell\) and any \(z,z_i\), \(\ell(z)\le\ell(z_i)+\ell'(z_i)(z-z_i)\). (Proof: by concavity, for \(0\lt\theta\le1\), \(\ell(z_i+\theta(z-z_i))\)\({}\ge(1-\theta)\ell(z_i)+\theta\ell(z)\); rearrange to \(\{\ell(z_i+\theta(z-z_i))-\ell(z_i)\}/\theta\)\({}\ge\ell(z)-\ell(z_i)\) and let \(\theta\to0\); the left side tends to \(\ell'(z_i)(z-z_i)\).) Taking the minimum over \(i\) keeps the inequality, and \(u(z_i)\le t_i(z_i)=\ell(z_i)\) gives equality at the grid points.

(b) Because \(\ell'\) is nonincreasing, the slopes satisfy \(\ell'(z_1)\ge\ell'(z_2)\ge\dots\). Two lines with slopes \(\ell'(z_i)\gt\ell'(z_{i+1})\) cross once; to the left of the crossing the line with the larger slope is lower. So \(t_i\le t_{i+1}\) exactly for \(z\le\widehat z_{i,i+1}\), and the crossing points are ordered, \(\widehat z_{i-1,i}\le\widehat z_{i,i+1}\) (a standard property of the lower envelope of lines with decreasing slopes). Hence the minimum is \(t_i\) between consecutive crossing points. The formula for \(\widehat z_{i,i+1}\) is the solution of \(t_i(z)=t_{i+1}(z)\), a linear equation in \(z\). (If \(\ell'(z_i)=\ell'(z_{i+1})\), then \(\ell\) is linear between the points and either tangent can be used.)

(c) On piece \(i\), \(e^{t_i(z)}\)\({}=\widetilde p(z_i)\exp\{-\lambda_i(z-z_i)\}\), which has the form (11.17) with \(k_i\lambda_i\) proportional to \(\widetilde p(z_i)\). Integrating each piece, \(\int_{a}^{b}e^{-\lambda(z-z_i)}\d z\)\({}=\lambda^{-1}[e^{-\lambda(a-z_i)}-e^{-\lambda(b-z_i)}]\), gives the pieces of \(C\). The two unbounded pieces are finite because \(\lambda_1=-\ell'(z_1)\lt0\) on \((-\infty,\widehat z_{1,2}]\) and \(\lambda_n\gt0\) on \((\widehat z_{n-1,n},\infty)\): this is why the outermost grid points must lie on either side of the mode. Dividing \(e^{u}\) by \(C\) gives \(q\) with \(k_i=\widetilde p(z_i)/(\lambda_iC)\). Since \(u\) is continuous (a minimum of continuous functions), so is \(q\): at \(\widehat z_{i,i+1}\), \(k_i\lambda_ie^{-\lambda_i(\widehat z_{i,i+1}-z_i)}\)\({}=k_{i+1}\lambda_{i+1}e^{-\lambda_{i+1}(\widehat z_{i,i+1}-z_{i+1})}\). These \(n-1\) continuity equations fix the ratios of the \(k_i\lambda_i\), and \(\int q=1\) fixes their common scale, which is the route of Exercise 11.8. Finally, \(e^{u}=Cq\ge\widetilde p\), so rejection sampling with proposal \(q\) and constant \(k=C\) is valid. ∎

Erratum. Some printings print the range in (11.17) and the following sentence incorrectly; the corrected form is the one above: the piece for \(z_i\) runs between \(\widehat z_{i-1,i}\) and \(\widehat z_{i,i+1}\), the intersections of the tangent at \(z_i\) with its neighbours.

Variants p. 531: a derivative-free version uses chords through pairs of grid points instead of tangents (Gilks, 1992), and adaptive rejection Metropolis sampling (Gilks, Best and Tan, 1995) handles densities that are not log-concave by following each rejection step with a Metropolis-Hastings step (Section 11.2.2).

Rejection sampling in high dimensions pp. 531–532

To see what dimension does to rejection sampling, the book samples \(p=\N(\mathbf 0,\sigma_p^2\mathbf I)\) in \(D\) dimensions using the proposal \(q=\N(\mathbf 0,\sigma_q^2\mathbf I)\), a nearly perfect proposal.

Proof

The ratio of the two densities is

\[ \frac{p(\mathbf z)}{q(\mathbf z)}=\frac{(2\pi\sigma_p^2)^{-D/2}e^{-\lVert\mathbf z\rVert^2/2\sigma_p^2}}{(2\pi\sigma_q^2)^{-D/2}e^{-\lVert\mathbf z\rVert^2/2\sigma_q^2}}=\Big(\frac{\sigma_q}{\sigma_p}\Big)^D\exp\Big\{-\frac{\lVert\mathbf z\rVert^2}{2}\Big(\frac1{\sigma_p^2}-\frac1{\sigma_q^2}\Big)\Big\}. \]

If \(\sigma_q\lt\sigma_p\), the exponent is positive and the ratio is unbounded as \(\lVert\mathbf z\rVert\to\infty\), so no \(k\) works. If \(\sigma_q\ge\sigma_p\), the exponent is \(\le0\) and the supremum, attained at \(\mathbf z=\mathbf 0\), is \((\sigma_q/\sigma_p)^D\). Both densities are normalised (\(Z_p=1\)), so by (11.14) the acceptance probability is \(1/k\). Numerically \(1000\ln1.01=9.950\), and \(e^{-9.950}=4.77\times10^{-5}=1/20\,959\), in line with the book's "approximately 1/20 000". The loss comes from the peak: in \(D\) dimensions, a 1% mismatch of scale is a factor \(1.01^D\) mismatch of height at the mode. ∎

The exponential loss is generic: an envelope must dominate \(\widetilde p\) everywhere, and in high dimensions any mismatch is multiplied across dimensions. Rejection sampling is useful in one or two dimensions and as a building block (for example inside Gibbs sampling) but not on its own in high dimensions p. 532.

11.1.4 Importance sampling

pp. 532–534 Importance sampling estimates expectations directly, without producing draws from \(p\). A grid approximation

\[ \E[f]\simeq\sum_{l=1}^{L}p(\mathbf z^{(l)})\,f(\mathbf z^{(l)}) \tag{11.18} \]

(times the volume of a grid cell, which the book absorbs) needs a number of points exponential in the dimension and wastes nearly all of them where \(p\) is negligible. Instead, draw \(\mathbf z^{(l)}\) from a proposal \(q\) and correct by the ratio \(p/q\):

\[ \E[f]=\int f(\mathbf z)p(\mathbf z)\,\d\mathbf z=\int f(\mathbf z)\frac{p(\mathbf z)}{q(\mathbf z)}q(\mathbf z)\,\d\mathbf z\simeq\frac1L\sum_{l=1}^{L}\frac{p(\mathbf z^{(l)})}{q(\mathbf z^{(l)})}f(\mathbf z^{(l)}). \tag{11.19} \]

The \(r_l=p(\mathbf z^{(l)})/q(\mathbf z^{(l)})\) are the importance weights. All samples are kept.

Proof

Apply Result 11.1 to the function \(\mathbf z\mapsto f(\mathbf z)p(\mathbf z)/q(\mathbf z)\) under \(q\). Its mean is \(\int f\frac pq\,q=\int fp=\E[f]\) (this is where positivity of \(q\) is needed: on a set where \(q=0\) but \(fp\ne0\), the integral of \(fp\) would be missed). Its second moment is \(\int f^2\frac{p^2}{q^2}q=\int f^2p^2/q\), which gives the variance.

For the optimum, by the Cauchy–Schwarz inequality, for any \(q\),

\[ \Big(\int|f|p\Big)^2=\Big(\int\frac{|f|p}{\sqrt q}\sqrt q\Big)^2\le\int\frac{f^2p^2}{q}\int q=\int\frac{f^2p^2}{q}, \]

with equality iff \(|f|p/\sqrt q\) is proportional to \(\sqrt q\), i.e. \(q\propto|f|p\). So the second moment is minimised by \(q^\star\), where it equals \((\E|f|)^2\); if \(f\ge0\), that is \(\E[f]^2\) and the variance is zero. The optimum needs \(\int|f|p\), the very kind of integral we are trying to compute, so it is a guide (put proposal mass where \(|f|p\) is large, as the book says on p. 532) rather than a recipe. ∎

Unnormalised densities p. 533. Usually \(p=\widetilde p/Z_p\) with \(Z_p\) unknown, and perhaps \(q=\widetilde q/Z_q\) too. Then

\[ \E[f]=\int f(\mathbf z)p(\mathbf z)\,\d\mathbf z=\frac{Z_q}{Z_p}\int f(\mathbf z)\frac{\widetilde p(\mathbf z)}{\widetilde q(\mathbf z)}q(\mathbf z)\,\d\mathbf z\simeq\frac{Z_q}{Z_p}\,\frac1L\sum_{l=1}^{L}\widetilde r_lf(\mathbf z^{(l)}), \tag{11.20} \]

with \(\widetilde r_l\)\({}=\widetilde p(\mathbf z^{(l)})/\widetilde q(\mathbf z^{(l)})\). The same samples estimate the unknown ratio:

\[ \frac{Z_p}{Z_q}=\frac1{Z_q}\int\widetilde p(\mathbf z)\,\d\mathbf z=\int\frac{\widetilde p(\mathbf z)}{\widetilde q(\mathbf z)}q(\mathbf z)\,\d\mathbf z\simeq\frac1L\sum_{l=1}^{L}\widetilde r_l, \tag{11.21} \]

and dividing (11.20) by (11.21) gives the self-normalised estimator

\[ \E[f]\simeq\sum_{l=1}^{L}w_lf(\mathbf z^{(l)}), \tag{11.22} \] \[ w_l=\frac{\widetilde r_l}{\sum_m\widetilde r_m}=\frac{\widetilde p(\mathbf z^{(l)})/q(\mathbf z^{(l)})}{\sum_m\widetilde p(\mathbf z^{(m)})/q(\mathbf z^{(m)})}. \tag{11.23} \]

Proof

(a) Write \(\widehat f_{\rm SN}=A_L/B_L\) with \(A_L\)\({}=\frac1L\sum_l\widetilde r_lf(\mathbf z^{(l)})\) and \(B_L=\frac1L\sum_l\widetilde r_l\), where here \(\widetilde r_l=\widetilde p/q\). By the strong law, \(A_L\to\int f\widetilde p=Z_p\mu\) and \(B_L\to\int\widetilde p=Z_p\gt0\) almost surely, so the ratio tends to \(\mu\). Replacing \(q\) by \(\widetilde q=Z_qq\) multiplies every \(\widetilde r_l\) by the same constant \(1/Z_q\), which cancels in (11.23).

(b, c) Linearise the ratio around the limits (the delta method). With \(a=Z_p\mu\), \(b=Z_p\),

\[ \frac{A_L}{B_L}-\mu=\frac{A_L-\mu B_L}{B_L}=\frac{1}{b}(A_L-\mu B_L)\Big\{1-\frac{B_L-b}{b}+O_p(L^{-1})\Big\}. \]

The leading term \((A_L-\mu B_L)/b\)\({}=\frac1L\sum_l\frac{\widetilde r_l}{Z_p}\{f(\mathbf z^{(l)})\)\({}-\mu\}\)\({}=\frac1L\sum_lw(\mathbf z^{(l)})\{f(\mathbf z^{(l)})\)\({}-\mu\}\) is an average of independent terms with mean \(\int q\,w(f-\mu)=\int p(f-\mu)=0\) and variance \(\E_q[w^2(f-\mu)^2]\), which gives (c). The next term, \(-(A_L-\mu B_L)(B_L-b)/b^2\), is a product of two averages each of order \(L^{-1/2}\); its expectation is \(-\cov(A_L-\mu B_L,B_L)/b^2\)\({}=-\frac{1}{L}\E_q[w^2(f-\mu)]\), of order \(1/L\), which is the bias in (b). Finally, \(\E_q[w^2g]=\int\frac{p^2}{q^2}g\,q=\E_p[wg]\). ∎

When importance sampling fails pp. 533–534

If \(p\,f\) puts much of its mass where \(q\) is small, a few weights dominate and the effective number of samples is far below \(L\). Worse, if no sample has landed there yet, the weights may look harmless and the estimate may still be badly wrong, "with no diagnostic indication" (p. 534). The cleanest way to see this is through the tails.

Proof

The approximation (Kong, 1992; Kong, Liu and Wong, 1994). By Result 11.11c, the variance is about \(\frac1L\E_p[w(f-\mu)^2]\), while independent sampling from \(p\) gives \(\frac1L\E_p[(f-\mu)^2]\). Write \(\E_p[w(f-\mu)^2]\)\({}=\E_p[w]\,\E_p[(f-\mu)^2]\)\({}+\cov_p\big(w,(f-\mu)^2\big)\). The rule of thumb drops the covariance, which is exact when \(f\) does not vary with the weight and reasonable when \(f\) is not concentrated where \(w\) is extreme; it is a property of \(q\) alone, which is its appeal (and its limit: it can be optimistic for a particular \(f\)). Then the inflation factor is \(\E_p[w]=\int p^2/q=\E_q[w^2]\). Since \(\E_q[w]=\int p=1\), \(\E_q[w^2]=1+\var_q[w]\ge1\), with equality iff \(w\) is constant, i.e. \(q=p\) (equivalently, by Jensen's inequality \(\E_q[w^2]\ge(\E_qw)^2\); toolkit).

The estimate. Replace \(\E_q[w^2]\)\({}=\E_q[\widetilde r^2]/(\E_q\widetilde r)^2\) (the normalising constants cancel) by sample averages: \(\frac{\frac1L\sum\widetilde r_l^2}{(\frac1L\sum\widetilde r_l)^2}\), so \(\widehat L_{\rm eff}\)\({}=L\big/\big(L\sum\widetilde r_l^2/(\sum\widetilde r_l)^2\big)\)\({}=(\sum\widetilde r_l)^2/\sum\widetilde r_l^2\)\({}=1/\sum w_l^2\). Bounds. By Cauchy–Schwarz, \((\sum w_l)^2\le L\sum w_l^2\), so \(1/\sum w_l^2\le L\) with equality iff all \(w_l=1/L\); and \(\sum w_l^2\le(\max_lw_l)\sum_lw_l\le1\) gives \(\widehat L_{\rm eff}\ge1\). ∎

Importance sampling on graphical models p. 534, (11.24)

With evidence \(\mathbf x\) on some nodes of a directed graph, the target is the posterior of the other nodes, \(p(\mathbf z\given\mathbf x)\)\({}\propto\widetilde p(\mathbf z)\), where \(\widetilde p\) is the joint (11.4) evaluated at the observed values. Uniform sampling draws the unobserved nodes uniformly over their states; then \(\widetilde q\) is constant and the weights are proportional to \(\widetilde p(\mathbf z)\). Likelihood weighted sampling (Fung and Chang, 1990; Shachter and Peot, 1990) runs ancestral sampling but clamps each evidence node at its observed value instead of sampling it. Its weight is

\[ r(\mathbf z)=\prod_{z_i\notin e}\frac{p(z_i\given\mathrm{pa}_i)}{p(z_i\given\mathrm{pa}_i)}\prod_{z_i\in e}\frac{p(z_i\given\mathrm{pa}_i)}{1}=\prod_{z_i\in e}p(z_i\given\mathrm{pa}_i). \tag{11.24} \]

11.1.5 Sampling-importance-resampling

pp. 534–536 Rejection sampling needs a constant \(k\), which may be unknown or so large that almost nothing is accepted. Sampling-importance-resampling (SIR; Rubin, 1987, 1988) avoids \(k\):

Proof

Given the first stage, the resampled value equals \(z^{(l)}\) with probability \(w_l\), so its distribution function is the sum of the weights of the points \(\le a\), which is (11.25) after inserting (11.23). Numerator and denominator divided by \(L\) are sample means of \(I(z\le a)\widetilde p/q\) and of \(\widetilde p/q\) under \(q\); both have finite means (\(\int I\,\widetilde p\le Z_p\) and \(Z_p\)), so by the strong law they converge almost surely to the integrals in (11.26), the denominator to \(Z_p\gt0\). The ratio converges to their ratio. Convergence of distribution functions at every \(a\) is convergence in distribution; unconditionally, \(\Pr(z^\ast\le a)\) is the expectation of (11.25), which converges too by bounded convergence (the ratio is between 0 and 1). For a fixed \(L\) the resampled values are only approximately from \(p\): their distribution is the self-normalised estimator of Result 11.11, with the same \(O(1/L)\) bias. If \(q=p\), all weights equal \(1/L\) and resampling returns draws from \(p\). ∎

For moments, resampling only adds noise: use the weighted original sample directly p. 536,

\[ \E[f(\mathbf z)]=\frac{\int f(\mathbf z)[\widetilde p(\mathbf z)/q(\mathbf z)]q(\mathbf z)\,\d\mathbf z}{\int[\widetilde p(\mathbf z)/q(\mathbf z)]q(\mathbf z)\,\d\mathbf z}\simeq\sum_{l=1}^{L}w_lf(\mathbf z^{(l)}). \tag{11.27} \]

(Given the first stage, the average of \(f\) over the resampled values has conditional mean \(\sum_lw_lf(\mathbf z^{(l)})\) and a positive conditional variance, so by the law of total variance it is never more precise than (11.27).) Resampling is useful when unweighted draws are needed, as in particle filters (Section 13.3.4).

Erratum. Some printings write \(\mathbf z_l\) in the last line of (11.27); it should be \(\mathbf z^{(l)}\).

11.1.6 Sampling and the EM algorithm

pp. 536–537 When the E step of EM (Chapter 9) cannot be done in closed form, it can be done by sampling. With hidden variables \(\mathbf Z\), data \(\mathbf X\) and parameters \(\bm\theta\), the M step maximises

\[ Q(\bm\theta,\bm\theta^{\rm old})=\int p(\mathbf Z\given\mathbf X,\bm\theta^{\rm old})\ln p(\mathbf Z,\mathbf X\given\bm\theta)\,\d\mathbf Z, \tag{11.28} \]

and Monte Carlo EM (Wei and Tanner, 1990) replaces it by an average over draws \(\mathbf Z^{(l)}\) from \(p(\mathbf Z\given\mathbf X,\bm\theta^{\rm old})\):

\[ Q(\bm\theta,\bm\theta^{\rm old})\simeq\frac1L\sum_{l=1}^{L}\ln p(\mathbf Z^{(l)},\mathbf X\given\bm\theta). \tag{11.29} \]

Data augmentation p. 537. For a full Bayesian treatment we want draws from \(p(\bm\theta\given\mathbf X)\) (or from the joint \(p(\bm\theta,\mathbf Z\given\mathbf X)\)), and suppose that the complete-data posterior \(p(\bm\theta\given\mathbf Z,\mathbf X)\) is easy to sample. The IP algorithm of Tanner and Wong (1987) alternates two steps, built on the identities

\[ p(\mathbf Z\given\mathbf X)=\int p(\mathbf Z\given\bm\theta,\mathbf X)\,p(\bm\theta\given\mathbf X)\,\d\bm\theta, \tag{11.30} \] \[ p(\bm\theta\given\mathbf X)=\int p(\bm\theta\given\mathbf Z,\mathbf X)\,p(\mathbf Z\given\mathbf X)\,\d\mathbf Z. \tag{11.31} \]

Proof

The map. The I step draws \(\bm\theta'\sim g\) and \(\mathbf Z\)\({}\sim p(\mathbf Z\given\bm\theta',\mathbf X)\), so each \(\mathbf Z^{(l)}\) has density \(\int p(\mathbf Z\given\bm\theta',\mathbf X)g(\bm\theta')\,\d\bm\theta'\), the right side of (11.30) with \(g\) in place of the posterior. As \(L\to\infty\), the mixture (11.32) tends to the expectation of \(p(\bm\theta\given\mathbf Z,\mathbf X)\) under this density (law of large numbers), which is \(\int K(\bm\theta',\bm\theta)g(\bm\theta')\d\bm\theta'\).

(a) If \(g=p(\cdot\given\mathbf X)\), the I-step density of \(\mathbf Z\) is \(p(\mathbf Z\given\mathbf X)\) by (11.30), and then the P step returns \(p(\bm\theta\given\mathbf X)\) by (11.31). The two identities are the sum rule applied to the joint \(p(\bm\theta,\mathbf Z\given\mathbf X)\), factorised by the product rule in the two possible orders.

(b) By Bayes' rule, \(p(\mathbf Z\given\bm\theta',\mathbf X)p(\bm\theta'\given\mathbf X)\)\({}=p(\bm\theta'\given\mathbf Z,\mathbf X)p(\mathbf Z\given\mathbf X)\). Hence

\[ p(\bm\theta'\given\mathbf X)K(\bm\theta',\bm\theta)=\int p(\bm\theta\given\mathbf Z,\mathbf X)\,p(\bm\theta'\given\mathbf Z,\mathbf X)\,p(\mathbf Z\given\mathbf X)\,\d\mathbf Z, \]

which is symmetric in \(\bm\theta\) and \(\bm\theta'\): this is detailed balance (11.40), and it implies (a) again (Result 11.16 below).

(c) With \(\pi=p(\cdot\given\mathbf X)\) and (a), \(g_{\rm new}-\pi\)\({}=\int K(\bm\theta',\cdot)\{g(\bm\theta')\)\({}-\pi(\bm\theta')\}\d\bm\theta'\). Taking absolute values inside and using \(\int K(\bm\theta',\bm\theta)\,\d\bm\theta=1\): \(\int|g_{\rm new}-\pi|\)\({}\le\int\!\int K(\bm\theta',\bm\theta)|g(\bm\theta')\)\({}-\pi(\bm\theta')|\,\d\bm\theta'\d\bm\theta\)\({}=\int|g-\pi|\). Tanner and Wong (1987) show that under mild positivity conditions the distance in fact tends to 0.

(d) With \(L=1\), the P step's "mixture" is the single density \(p(\bm\theta\given\mathbf Z^{(1)},\mathbf X)\), and the next I step draws \(\bm\theta\) from it and then \(\mathbf Z\) given \(\bm\theta\): alternately drawing \(\bm\theta\given\mathbf Z\) and \(\mathbf Z\given\bm\theta\) from their full conditionals is Gibbs sampling. ∎

The distinction between parameters and hidden variables is artificial here p. 537; from now on the problem is simply to sample a given distribution \(p(\mathbf z)\).

11.2 Markov chain Monte Carlo

pp. 537–539 Rejection and importance sampling degrade exponentially with dimension because a single fixed proposal must cover the whole of \(p\). Markov chain Monte Carlo (MCMC) instead lets the proposal depend on the current state: from \(\mathbf z^{(\tau)}\) it proposes a candidate \(\mathbf z^\star\)\({}\sim q(\mathbf z\given\mathbf z^{(\tau)})\), a local move, and accepts or rejects it. The successive states form a Markov chain whose distribution, under conditions made precise in Section 11.2.1, converges to \(p\). As before, only \(\widetilde p=Z_pp\) needs to be evaluated. The method comes from physics (Metropolis and Ulam, 1949; Metropolis et al., 1953) and entered statistics around 1990.

In the Metropolis algorithm the proposal is symmetric, \(q(\mathbf z_A\given\mathbf z_B)\)\({}=q(\mathbf z_B\given\mathbf z_A)\) for all \(\mathbf z_A,\mathbf z_B\), and the candidate is accepted with probability

\[ A(\mathbf z^\star,\mathbf z^{(\tau)})=\min\Big(1,\frac{\widetilde p(\mathbf z^\star)}{\widetilde p(\mathbf z^{(\tau)})}\Big). \tag{11.33} \]

Repeating the current state after a rejection is essential: dropping repeats would over-represent the states that are easy to leave. The draws are dependent, so Result 11.2 applies; the book suggests keeping every \(M\)th draw, which gives nearly independent draws but, as noted there, does not make estimates from a given run more precise.

Random walks p. 539. To see why local moves are slow, take the random walk on the integers with

\[ p(z^{(\tau+1)}=z^{(\tau)})=0.5, \tag{11.34} \] \[ p(z^{(\tau+1)}=z^{(\tau)}+1)=0.25, \tag{11.35} \] \[ p(z^{(\tau+1)}=z^{(\tau)}-1)=0.25, \tag{11.36} \]

started at \(z^{(0)}=0\).

Erratum. Some printings say \(z^{(1)}=0\) after (11.36); the walk starts at \(z^{(0)}=0\).

Proof

Write \(z^{(\tau)}=z^{(\tau-1)}+\xi_\tau\), where the step \(\xi_\tau\in\{-1,0,1\}\) has probabilities \(0.25,0.5,0.25\) and is independent of the past. Then \(\E[\xi_\tau]=0\) and \(\E[\xi_\tau^2]=0.25+0.25=\tfrac12\). By induction \(\E[z^{(\tau)}]=\E[z^{(\tau-1)}]=\dots=0\). Squaring, \(\E[(z^{(\tau)})^2]\)\({}=\E[(z^{(\tau-1)})^2]\)\({}+2\E[z^{(\tau-1)}]\E[\xi_\tau]+\E[\xi_\tau^2]\)\({}=\E[(z^{(\tau-1)})^2]+\tfrac12\), using independence of \(\xi_\tau\) from \(z^{(\tau-1)}\). With \(\E[(z^{(0)})^2]=0\), induction gives \(\tau/2\). (In my simulation of 200 000 walks, the mean of \((z^{(50)})^2\) was 25.03.) By the central limit theorem, \(z^{(\tau)}/\sqrt{\tau/2}\) is approximately standard normal: to travel a distance \(d\), a random walk needs of order \(d^2\) steps. ∎

This square-root law is the central inefficiency of MCMC with local proposals, and avoiding it motivates over-relaxation (Section 11.3) and hybrid Monte Carlo (Section 11.5).

11.2.1 Markov chains

pp. 539–541 We now make precise when a Markov chain "converges to \(p\)".

In words, (11.38) says that the probability of being at \(\mathbf z\) next is the probability of being at each \(\mathbf z'\) now times the probability of moving from there to \(\mathbf z\), summed over \(\mathbf z'\). Detailed balance says more: in equilibrium the probability flow from \(\mathbf z\) to \(\mathbf z'\) equals the flow back, for each pair separately. Equivalently, if \(\mathbf z^{(m)}\sim p^\star\), the pair \((\mathbf z^{(m)},\mathbf z^{(m+1)})\) has the same joint distribution as \((\mathbf z^{(m+1)},\mathbf z^{(m)})\): a film of the stationary chain looks the same run backwards, hence "reversible".

Proof

The first equality replaces each term by the other side of (11.40); the second takes \(p^\star(\mathbf z)\) out of the sum; the third uses that the transition probabilities out of \(\mathbf z\) sum to one. In the continuous case, integrate both sides of (11.40) over \(\mathbf z'\) in the same way.

The converse fails. On three states, let the chain move deterministically \(1\to2\to3\to1\). The uniform distribution is invariant (each state receives exactly the mass of its predecessor), but \(p^\star(1)T(1,2)=\tfrac13\ne0=p^\star(2)T(2,1)\). A chain can circulate in equilibrium; detailed balance forbids net circulation. So detailed balance is sufficient but not necessary, as the book says. ∎

Ergodicity p. 540. Invariance is not enough: the identity transition leaves every distribution invariant and never moves. We also need the chain to forget its starting point: \(p(\mathbf z^{(m)})\to p^\star\) as \(m\to\infty\) for every initial distribution. Such a chain is ergodic, and \(p^\star\) is then its unique equilibrium distribution. Even detailed balance plus "can reach every state" is not enough: the two-state chain that always switches state satisfies detailed balance with respect to \((\tfrac12,\tfrac12)\), but started in state 1 it is in state 1 at every even step, so \(p(\mathbf z^{(m)})\) oscillates and never converges (periodicity). The book quotes Neal (1993) for "weak restrictions" that suffice; for finite state spaces a complete proof is short.

Proof

Treat distributions as row vectors, so one step maps \(\mu\) to \(\mu T\) (this is (11.38)), and write \(\lVert d\rVert_1=\sum|d(\mathbf z)|\).

Step 1 (contraction). Let \(P=T^m\) and \(c(\mathbf z')\)\({}=\min_{\mathbf z}P(\mathbf z,\mathbf z')\), so \(\beta=\sum_{\mathbf z'}c(\mathbf z')\). For two distributions \(\mu,\nu\), let \(d=\mu-\nu\); then \(\sum_{\mathbf z}d(\mathbf z)=0\), so subtracting \(c(\mathbf z')\sum_{\mathbf z}d(\mathbf z)=0\) changes nothing:

\[ (dP)(\mathbf z')=\sum_{\mathbf z}d(\mathbf z)\big\{P(\mathbf z,\mathbf z')-c(\mathbf z')\big\}. \]

The braces are \(\ge0\) and, summed over \(\mathbf z'\), equal \(1-\beta\). Hence

\[ \lVert dP\rVert_1\le\sum_{\mathbf z'}\sum_{\mathbf z}|d(\mathbf z)|\big\{P(\mathbf z,\mathbf z')-c(\mathbf z')\big\}=(1-\beta)\lVert d\rVert_1 . \]

The same argument with \(c\equiv0\) shows \(\lVert dT\rVert_1\le\lVert d\rVert_1\): one step of any chain never increases the distance between two distributions.

Step 2 (existence and uniqueness). The set of distributions is closed and bounded in \(\R^{|S|}\), hence complete, and \(\mu\mapsto\mu P\) is a contraction on it (factor \(1-\beta\lt1\)). By the Banach fixed-point theorem it has a unique fixed point \(\pi\). Then \(\pi T\) is also fixed by \(P\), since \((\pi T)P=(\pi P)T=\pi T\); by uniqueness \(\pi T=\pi\), so \(\pi\) is invariant for \(T\). Any invariant distribution of \(T\) is invariant for \(P\), so equals \(\pi\). Call it \(p^\star\). (In MCMC we know \(p^\star\) in advance; the point of (a) is that there is no other.)

Step 3 (convergence). Write \(n=km+j\) with \(0\le j\lt m\). Since \(p^\star=p^\star T^n\), \(\lVert\mu T^n-p^\star\rVert_1\)\({}=\lVert(\mu-p^\star)P^kT^j\rVert_1\)\({}\le(1-\beta)^k\lVert\mu-p^\star\rVert_1\)\({}\le2(1-\beta)^k\), by Step 1 applied \(k\) times and the non-expansion for the remaining \(j\) steps.

Step 4 (ergodic theorem, sketch). Fix a state \(\mathbf a\). By the Markov property the chain restarts afresh at each visit to \(\mathbf a\), so the excursions between successive visits are independent and identically distributed; they have finite mean length because the geometric convergence of Step 3 makes long excursions exponentially unlikely. Applying the strong law of large numbers to the sums of \(f\) over excursions and to their lengths gives the limit \(\E[\text{sum of }f\text{ over an excursion}]/\E[\text{length}]\), and the same identity applied to the indicator of each state shows that this ratio is \(\E_{p^\star}[f]\). (Full details: Levin, Peres and Wilmer, 2017, Chapter 4.)

Positive power. For a finite chain, "some \(T^m\gt0\) entrywise" is equivalent to irreducibility (every state can reach every other) plus aperiodicity (the return times to a state have greatest common divisor 1); the number-theoretic argument is in Levin, Peres and Wilmer (2017), Proposition 1.7. ∎

General state spaces. For continuous \(\mathbf z\) the corresponding theorem (Tierney, 1994, Theorem 1) reads: if \(T\) has invariant distribution \(p^\star\), is \(p^\star\)-irreducible (from any start, every set of positive \(p^\star\)-probability is reached with positive probability) and aperiodic, then \(p(\mathbf z^{(m)})\to p^\star\) in total variation for \(p^\star\)-almost every starting point, and ergodic averages converge to \(\E_{p^\star}[f]\) for every \(f\) with \(\E_{p^\star}|f|\lt\infty\). These are the "weak restrictions" of the book. In practice: build \(T\) with \(p\) invariant (usually through detailed balance), and check that the chain can get from anywhere to anywhere.

Combining transitions pp. 540–541. Samplers are usually assembled from base transitions \(B_1,\dots,B_K\), either as a mixture (pick \(B_k\) with probability \(\alpha_k\))

\[ T(\mathbf z',\mathbf z)=\sum_{k=1}^{K}\alpha_kB_k(\mathbf z',\mathbf z),\qquad\alpha_k\ge0,\ \sum_k\alpha_k=1, \tag{11.42} \]

or as a composition (apply \(B_1\), then \(B_2\), ..., then \(B_K\))

\[ T(\mathbf z',\mathbf z)=\sum_{\mathbf z_1}\dots\sum_{\mathbf z_{K-1}}B_1(\mathbf z',\mathbf z_1)\dots B_{K-1}(\mathbf z_{K-2},\mathbf z_{K-1})B_K(\mathbf z_{K-1},\mathbf z). \tag{11.43} \]

Erratum. In some printings the last sum in (11.43) runs over \(\mathbf z_{n-1}\); it should be \(\mathbf z_{K-1}\).

Proof

Use matrix notation on a finite space (the general case is identical with integrals): \(p^\star\) is a row vector, invariance reads \(p^\star B=p^\star\), and with \(\mathbf D=\diag(p^\star)\) detailed balance reads \(\mathbf DB=B^\T\mathbf D\) (entry \((\mathbf z,\mathbf z')\) of both sides is (11.40)).

(a) Mixture: \(p^\star T\)\({}=\sum_k\alpha_kp^\star B_k\)\({}=\sum_k\alpha_kp^\star\)\({}=p^\star\). Composition: (11.43) is the matrix product \(T=B_1B_2\cdots B_K\), and \(p^\star B_1B_2\cdots B_K\)\({}=p^\star B_2\cdots B_K\)\({}=\dots\)\({}=p^\star\).

(b) \(\mathbf DT\)\({}=\sum_k\alpha_k\mathbf DB_k\)\({}=\sum_k\alpha_kB_k^\T\mathbf D\)\({}=T^\T\mathbf D\).

(c) For \(T=B_1B_2\): \(\mathbf DB_1B_2\)\({}=B_1^\T\mathbf DB_2\)\({}=B_1^\T B_2^\T\mathbf D\)\({}=(B_2B_1)^\T\mathbf D\). This equals \(T^\T\mathbf D=(B_1B_2)^\T\mathbf D\) only if \(B_1B_2=B_2B_1\), which fails in general. Example: the systematic-scan Gibbs sampler of Section 11.3 on two binary variables with a non-product distribution.

(d) For \(S=B_1\cdots B_KB_K\cdots B_1\), move \(\mathbf D\) from left to right through all \(2K\) factors, each time using \(\mathbf DB_k=B_k^\T\mathbf D\):

\[ \mathbf DS=B_1^\T\cdots B_K^\T B_K^\T\cdots B_1^\T\,\mathbf D=(B_1\cdots B_KB_K\cdots B_1)^\T\mathbf D=S^\T\mathbf D, \]

because the transpose of a product is the product of the transposes in reverse order, and the palindromic order is its own reverse. ∎

The typical use: each \(B_k\) changes only a subset of the variables (a "coordinate" move). A single \(B_k\) is then not ergodic, but the mixture or composition can be.

11.2.2 The Metropolis-Hastings algorithm

pp. 541–542 Hastings (1970) allowed asymmetric proposals. From the current state \(\mathbf z^{(\tau)}\), draw \(\mathbf z^\star\)\({}\sim q_k(\mathbf z\given\mathbf z^{(\tau)})\) and accept with probability

\[ A_k(\mathbf z^\star,\mathbf z^{(\tau)})=\min\Big(1,\frac{\widetilde p(\mathbf z^\star)\,q_k(\mathbf z^{(\tau)}\given\mathbf z^\star)}{\widetilde p(\mathbf z^{(\tau)})\,q_k(\mathbf z^\star\given\mathbf z^{(\tau)})}\Big), \tag{11.44} \]

where \(k\) labels one of several available transition types (to be mixed or composed as in (11.42)–(11.43)). The unknown \(Z_p\) cancels, and for a symmetric proposal the \(q\) factors cancel too, giving (11.33). The extra factor corrects for a proposal that goes more easily one way than the other.

Proof

The transition. From \(\mathbf z\), the chain moves to \(\mathbf z'\ne\mathbf z\) with density \(q_k(\mathbf z'\given\mathbf z)A_k(\mathbf z',\mathbf z)\) (propose it, then accept it), and stays at \(\mathbf z\) with the remaining probability \(r(\mathbf z)\)\({}=1\)\({}-\int q_k(\mathbf z'\given\mathbf z)A_k(\mathbf z',\mathbf z)\,\d\mathbf z'\). So

\[ T(\mathbf z,\d\mathbf z')=q_k(\mathbf z'\given\mathbf z)A_k(\mathbf z',\mathbf z)\,\d\mathbf z'+r(\mathbf z)\,\delta_{\mathbf z}(\d\mathbf z'). \]

Off-diagonal part. Multiply \(A_k(\mathbf z',\mathbf z)=\min(1,R)\), with \(R\)\({}=\dfrac{p(\mathbf z')q_k(\mathbf z\given\mathbf z')}{p(\mathbf z)q_k(\mathbf z'\given\mathbf z)}\), by \(p(\mathbf z)q_k(\mathbf z'\given\mathbf z)\gt0\): since \(c\min(1,R)=\min(c,cR)\) for \(c\gt0\), the first line of (11.45) follows. (\(R\) uses \(p\), but \(p/\widetilde p=1/Z_p\) cancels, so (11.44) is the same quantity.) The second line just swaps the arguments of the symmetric function \(\min\), and the third reads the first line backwards with \(\mathbf z\) and \(\mathbf z'\) exchanged. If \(p(\mathbf z)q_k(\mathbf z'\given\mathbf z)=0\), the move \(\mathbf z\to\mathbf z'\) never happens from a state of positive probability, and both sides of (11.45) are 0 with the convention \(A=0\) when the numerator vanishes (a move to a state of zero probability is rejected).

Diagonal part. The point mass contributes \(p(\mathbf z)r(\mathbf z)\) on the diagonal \(\mathbf z'=\mathbf z\) to both \(p(\d\mathbf z)T(\mathbf z,\d\mathbf z')\) and its mirror image, so it is symmetric automatically. Detailed balance holds, and invariance follows from Result 11.16.

Metropolis. If \(q_k(\mathbf a\given\mathbf b)\)\({}=q_k(\mathbf b\given\mathbf a)\), then \(R\)\({}=\widetilde p(\mathbf z')/\widetilde p(\mathbf z)\), which is (11.33).

Ergodicity. If \(q_k(\mathbf z'\given\mathbf z)\gt0\) for all \(\mathbf z,\mathbf z'\) (the book's sufficient condition, p. 538), then from any \(\mathbf z\) every set \(B\) with \(p(B)\gt0\) has \(T(\mathbf z,B)\ge\int_Bq_kA_k\gt0\) in one step, which gives irreducibility and aperiodicity, and Tierney's theorem applies. ∎

Erratum. Some printings have a different (incorrect) form of (11.45) and write \(\mathbf z^\tau\) for \(\mathbf z^{(\tau)}\) just before (11.44); the corrected forms are given above.

Other acceptance functions also give detailed balance (for example Barker's \(R/(1+R)\)), but \(\min(1,R)\) accepts as often as possible and, by Peskun's (1973) theorem, gives estimates with the smallest asymptotic variance among them for the same proposal.

Choosing the proposal scale pp. 541–542

With a Gaussian random-walk proposal \(\mathbf z^\star\)\({}=\mathbf z^{(\tau)}+\rho\,\bm\epsilon\), \(\bm\epsilon\sim\N(\mathbf 0,\mathbf I)\), a small \(\rho\) accepts nearly everything but moves slowly; a large \(\rho\) proposes into regions of low probability and is mostly rejected. On a Gaussian with very different standard deviations \(\sigma_{\max}\) and \(\sigma_{\min}\) along its principal axes (Figure 11.10), the scale must be of order \(\sigma_{\min}\) to keep the acceptance rate reasonable, and the chain must then random-walk along the long axis.

11.3 Gibbs sampling

pp. 542–546 Gibbs sampling (Geman and Geman, 1984) updates one variable at a time by drawing it from its full conditional \(p(z_i\given\mathbf z_{\setminus i})\), the distribution of \(z_i\) given all the other variables at their current values. For three variables, one sweep draws \(z_1^{(\tau+1)}\) from

\[ p(z_1\given z_2^{(\tau)},z_3^{(\tau)}), \tag{11.46} \] \[ \text{then }z_2^{(\tau+1)}\text{ from }p(z_2\given z_1^{(\tau+1)},z_3^{(\tau)}), \tag{11.47} \] \[ \text{then }z_3^{(\tau+1)}\text{ from }p(z_3\given z_1^{(\tau+1)},z_2^{(\tau+1)}), \tag{11.48} \]

always using the newest values.

Erratum. Some printings write \(p(z_i\given\{\mathbf z_{\setminus i})\) on p. 544 (a stray brace); it is \(p(z_i\given\mathbf z_{\setminus i})\).

Proof

(a) Use \(p(\mathbf z)\)\({}=p(z_k\given\mathbf z_{\setminus k})p(\mathbf z_{\setminus k})\) for both \(\mathbf z\) and \(\mathbf z^\star\), and the proposal densities as defined. Because \(\mathbf z^\star_{\setminus k}\)\({}=\mathbf z_{\setminus k}\), every factor in the numerator of (11.49) matches one in the denominator: \(p(z_k^\star\given\mathbf z^\star_{\setminus k})\) with \(p(z_k^\star\given\mathbf z_{\setminus k})\), \(p(\mathbf z^\star_{\setminus k})\) with \(p(\mathbf z_{\setminus k})\), and \(p(z_k\given\mathbf z^\star_{\setminus k})\) with \(p(z_k\given\mathbf z_{\setminus k})\). So \(A=\min(1,1)=1\): the proposal is always accepted, and the MH transition is exactly the Gibbs update.

(b) Directly (the solution of Exercise 11.11): \(B_k(\mathbf z,\mathbf z')\)\({}=p(z_k'\given\mathbf z_{\setminus k})\) if \(\mathbf z'_{\setminus k}=\mathbf z_{\setminus k}\), and 0 otherwise. For such a pair,

\[ p(\mathbf z)B_k(\mathbf z,\mathbf z')=p(z_k\given\mathbf z_{\setminus k})\,p(\mathbf z_{\setminus k})\,p(z'_k\given\mathbf z_{\setminus k}), \]

which is symmetric under exchanging \(z_k\) and \(z_k'\) (the common \(\mathbf z_{\setminus k}=\mathbf z'_{\setminus k}\) is untouched), so it equals \(p(\mathbf z')B_k(\mathbf z',\mathbf z)\). For pairs that differ outside coordinate \(k\), both sides are 0. Invariance then follows from Result 11.16. The book's own argument is the same fact read differently: the update leaves the marginal of \(\mathbf z_{\setminus k}\) untouched and redraws \(z_k\) from the correct conditional, so if \(\mathbf z\sim p\) before, the pair (marginal, conditional) and hence the joint is \(p\) after.

(c) By Result 11.18: invariance for both; detailed balance for the mixture (random scan) and for the palindromic composition \(B_1\cdots B_MB_M\cdots B_1\); not in general for \(B_1\cdots B_M\). ∎

Ergodicity p. 544. The book's sufficient condition is that no conditional distribution is zero anywhere. For a density on \(\R^M\), a convenient version is: if \(p(\mathbf z)\gt0\) for all \(\mathbf z\), then after one sweep the chain has moved from \(\mathbf z\) to \(\mathbf z'\) with density \(\prod_{j}p(z_j'\given z'_1,\dots,z'_{j-1},z_{j+1},\dots,z_M)\)\({}\gt0\) for every \(\mathbf z'\), so every set of positive probability is reached in one sweep from any starting point: the sweep is irreducible and aperiodic, and Tierney's theorem (Section 11.2.1) applies. When the support has a shape that the coordinate moves cannot cross, ergodicity can fail (Exercise 11.12, Figure 11.15: two regions that share no value of \(z_1\) and no value of \(z_2\), so the sampler can never leave the region it starts in). As with every MCMC method, the first part of the run (burn-in) is discarded, and the draws are correlated.

Gibbs sampling a correlated Gaussian pp. 544–545, Figure 11.11

The book argues that the step size of Gibbs sampling is set by the conditional widths \(l\), while the distribution extends a distance \(L\) along its long axis, so that a random walk needs about \((L/l)^2\) steps. For a bivariate Gaussian everything can be computed exactly.

Proof

(a) The conditional of one block of a Gaussian given the other (Section 2.3.1; toolkit) with covariance entries \(\Sigma_{11}=\sigma_1^2\), \(\Sigma_{12}=\rho\sigma_1\sigma_2\), \(\Sigma_{22}=\sigma_2^2\): the mean is \(\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(z_2-\mu_2)\) and the variance \(\Sigma_{11}-\Sigma_{12}^2/\Sigma_{22}\)\({}=\sigma_1^2(1-\rho^2)\).

(b) Standardise (\(\mu_i=0\), \(\sigma_i=1\)); Gibbs sampling commutes with this affine change of each coordinate. One sweep is \(z_1'=\rho z_2+s\nu_1\), \(z_2'=\rho z_1'+s\nu_2\), with \(s=\sqrt{1-\rho^2}\) and independent standard normals \(\nu_1,\nu_2\). Substituting the previous sweep's \(z_2=\rho z_1+s\nu_2^{\rm old}\): \(z_1'=\rho^2z_1+\rho s\nu^{\rm old}_2+s\nu_1\), an autoregression with coefficient \(\phi=\rho^2\) and noise independent of \(z_1\). In stationarity the lag-\(k\) autocorrelation is \(\phi^k\), and Result 11.2 gives \(\tau_{\rm int}=(1+\phi)/(1-\phi)\).

(c) The conditional variance is \(1-\rho^2\). The covariance matrix \(\begin{pmatrix}1&\rho\\\rho&1\end{pmatrix}\) has eigenvalues \(1\pm\rho\) with eigenvectors along the diagonals, so the long-axis standard deviation is \(\sqrt{1+\rho}\). Then \((L/l)^2=(1+\rho)/(1-\rho^2)=1/(1-\rho)\), and \(\tau_{\rm int}=(1+\rho^2)/\{(1-\rho)(1+\rho)\}\), whose ratio to \(1/(1-\rho)\) is \((1+\rho^2)/(1+\rho)\to1\). ∎

If the variables were uncorrelated, Gibbs sampling would produce independent draws at every sweep. For this toy problem a rotation would decorrelate them, but in real problems such a transformation is not available p. 544.

Over-relaxation p. 545, (11.50)

Adler's (1981) over-relaxation applies when the full conditionals are Gaussian. This is more general than a Gaussian joint: \(p(z,y)\propto\exp(-z^2y^2-z^2-y^2)\) is not Gaussian, but each conditional is (for example \(z\given y\sim\N(0,\{2(1+y^2)\}^{-1})\)). (The book's example \(\exp(-z^2y^2)\) has Gaussian conditionals but cannot be normalised, since \(\int\!\int e^{-z^2y^2}\d z\,\d y\)\({}=\int\sqrt\pi/|y|\,\d y\)\({}=\infty\); the extra terms fix this.) If the conditional of \(z_i\) has mean \(\mu_i\) and variance \(\sigma_i^2\), replace \(z_i\) by

\[ z_i'=\mu_i+\alpha(z_i-\mu_i)+\sigma_i(1-\alpha^2)^{1/2}\nu, \tag{11.50} \]

with \(\nu\sim\N(0,1)\) and \(-1\lt\alpha\lt1\). For \(\alpha=0\) this is the Gibbs update; for \(\alpha\lt0\) the new value is pushed to the other side of the conditional mean.

Erratum. Some printings write \(\alpha_i^2\) in the last term of (11.50); it is \(\alpha^2\).

Proof

(a) Given \(\mathbf z_{\setminus i}\), \(\mu_i\) and \(\sigma_i\) are constants and \(z_i'\) is a linear combination of the independent Gaussians \(z_i\) and \(\nu\), so it is Gaussian, with mean \(\mu_i+\alpha\cdot0+0=\mu_i\) and variance \(\alpha^2\sigma_i^2+\sigma_i^2(1-\alpha^2)\)\({}=\sigma_i^2\). (Only the first two moments of \(\nu\) are needed for the mean and variance, as in the exercise; Gaussianity of \(\nu\) gives the full distribution.) Since \(\mathbf z_{\setminus i}\) is unchanged and \(z_i'\) has the correct conditional distribution, the joint distribution after the update is again \(p\).

(b) Given \(\mathbf z_{\setminus i}\), the pair \((z_i,z_i')\) is bivariate Gaussian with both marginals \(\N(\mu_i,\sigma_i^2)\) and covariance \(\E[(z_i-\mu_i)(z_i'-\mu_i)]=\alpha\sigma_i^2\). A bivariate Gaussian with equal marginals is symmetric under exchanging its two coordinates, so \(p(z_i\given\mathbf z_{\setminus i})T(z_i,z_i')\)\({}=p(z_i'\given\mathbf z_{\setminus i})T(z_i',z_i)\); multiplying by \(p(\mathbf z_{\setminus i})\) gives (11.40). Hence the over-relaxed sweep is built from reversible pieces, as in Result 11.20.

(c) With standardised variables, each over-relaxed update is linear plus Gaussian noise, so a sweep is a linear autoregression \(\mathbf z'=\mathbf A\mathbf z+\bm\eta\). I computed its stationary covariance \(\mathbf S\) (from \(\mathbf S\)\({}=\mathbf A\mathbf S\mathbf A^\T+\cov[\bm\eta]\)) and the autocorrelations \(\mathbf e_1^\T\mathbf A^k\mathbf S\mathbf e_1/S_{11}\) exactly, then summed them as in Result 11.2; the \(\alpha=0\) value agrees with \((1+\rho^2)/(1-\rho^2)\). Values below 1 mean negatively correlated successive draws, which beat independent sampling for the mean; the gain is largest for linear functions and smaller for others such as variances. ∎

With \(\alpha\) close to \(-1\), each update reflects the variable to the other side of its conditional mean with little noise, so successive sweeps keep moving in the same direction along the long axis instead of diffusing: directed motion replaces the random walk. Ordered over-relaxation (Neal, 1999) achieves a similar effect for non-Gaussian conditionals: draw several values from the conditional, sort them together with the current value, and move to the value in the mirror-image rank.

Gibbs sampling on graphical models, ICM and blocking pp. 545–546

On a graphical model, the full conditional of a node involves only its Markov blanket: its neighbours in an undirected graph, and its parents, children and co-parents in a directed graph (the book's Figure 11.12).

Iterated conditional modes (ICM, Section 8.3.3) replaces each draw by the mode of the conditional: a greedy, deterministic relative of Gibbs sampling that climbs to a local maximum of \(p\) instead of sampling it. Blocking Gibbs (Jensen, Kong and Kjærulff, 1995) updates groups of variables jointly from their joint conditional; the proof of Result 11.20 goes through unchanged with \(z_k\) replaced by a block \(\mathbf z_B\) (marginal of the rest unchanged, block drawn from the correct conditional). Strongly dependent variables put in the same block no longer slow each other down; in the bivariate Gaussian, one block of both variables gives independent draws.

11.4 Slice sampling

pp. 546–548 Metropolis is sensitive to its step size. Slice sampling (Neal, 2003) adapts the step size automatically. It adds an auxiliary variable \(u\) and samples uniformly from the region under the graph of \(\widetilde p\):

\[ \widehat p(z,u)=\begin{cases}1/Z_p&\text{if }0\le u\le\widetilde p(z),\\ 0&\text{otherwise,}\end{cases} \tag{11.51} \]

with \(Z_p=\int\widetilde p(z)\,\d z\). Its marginal is

\[ \int\widehat p(z,u)\,\d u=\int_0^{\widetilde p(z)}\frac1{Z_p}\,\d u=\frac{\widetilde p(z)}{Z_p}=p(z), \tag{11.52} \]

so sampling \((z,u)\) and discarding \(u\) samples \(p\) (this is Result 11.6 again). Under \(\widehat p\), the conditional of \(u\) given \(z\) is uniform on \([0,\widetilde p(z)]\), and the conditional of \(z\) given \(u\) is uniform on the slice \(S_u=\{z:\widetilde p(z)\gt u\}\). Gibbs sampling of \((z,u)\) alternates the two (Figure 11.13a). The second draw is the hard part, since the slice may be a union of intervals whose ends are unknown; it is enough to replace it by any update of \(z\) that leaves the uniform distribution on \(S_u\) invariant.

Proof

Stepping out produces the same interval from \(z\) and from \(z'\). The random initial placement fixes a grid \(G=\{z_{\min}^{0}+jw:j\in\mathbb Z\}\) whose offset modulo \(w\) is uniform, whatever the current point: the initial interval is the grid cell containing the current point. Stepping out then moves the left end down the grid until it hits a grid point outside \(S\), and the right end up the grid likewise. So the final interval is \((L,R]\) where \(L\) is the largest grid point \(\le z\) outside \(S\) and \(R\) the smallest grid point \(\gt z\) outside \(S\), and every grid point strictly between \(L\) and \(R\) lies in \(S\). Now take any \(z'\in S\cap(L,R)\). Starting from \(z'\) with the same grid, the largest grid point \(\le z'\) outside \(S\) is again \(L\) (all grid points in \((L,z']\) are in \(S\)), and similarly for \(R\). Since the grid offset has the same distribution from either starting point, \(\Pr(\text{interval }I\given z)\)\({}=\Pr(\text{interval }I\given z')\) for every interval \(I\) containing both. (Without the random placement this fails: a fixed interval centred at \(z\) would make the interval depend on \(z\) in a non-symmetric way. The book's description omits this point.)

Shrinkage is symmetric. Fix the interval \(I\) and two points \(z,z'\in S\cap I\). A run from \(z\) that ends at \(z'\) consists of rejected points \(y_1,\dots,y_m\notin S\), each drawn uniformly from the current interval, followed by \(z'\). At each rejection the interval keeps the side containing the current point \(z\). Because the run ends by drawing \(z'\), every interval along the way contains \(z'\) too, so no \(y_j\) lies between \(z\) and \(z'\); hence "\(y_j\lt z\)" and "\(y_j\lt z'\)" agree, and starting from \(z'\) the same rejected points shrink the interval in exactly the same way, ending with an interval from which \(z\) can be drawn. The density of the run is the product of \(1/(\text{length of the current interval})\) over its draws, and these lengths are the same for the run from \(z\) and its mirror run from \(z'\). Summing (integrating) over all sequences of rejected points gives \(\Pr(z\to z'\given I)=\Pr(z'\to z\given I)\) as densities.

Conclusion. \(T_u(z,z')\)\({}=\sum_I\Pr(I\given z)\Pr(z\to z'\given I)\) is symmetric, so \(\frac{1}{|S|}T_u(z,z')=\frac1{|S|}T_u(z',z)\): detailed balance for the uniform distribution on \(S\), whose invariance follows (Result 11.16). Under \(\widehat p\), the conditional of \(z\) given \(u\) is exactly that uniform distribution, so the \(z\) update leaves \(\widehat p\) invariant; step 1 is an exact Gibbs draw of \(u\given z\), which also leaves it invariant. The composition leaves \(\widehat p\) invariant, and its \(z\)-marginal \(p\) (11.52). Shrinkage terminates with probability one when \(\widetilde p\) is continuous, because the interval always contains the current point, which lies in the open set \(S\). In my run on a two-component Gaussian mixture (200 000 updates), the draws matched the target's distribution function (Kolmogorov–Smirnov \(p=0.13\) on every 50th draw). ∎

The width \(w\) affects only efficiency: if it is too small, stepping out takes more steps (linear in the slice width divided by \(w\)); if too large, shrinkage takes more rejections (logarithmic in \(w\) divided by the slice width). Neal (2003) also describes a doubling procedure and limits on the number of steps, which need an extra acceptance check.

11.5 The hybrid Monte Carlo algorithm

p. 548 Random-walk proposals travel a distance proportional to the square root of the number of steps, and larger steps are rejected more often. Hybrid Monte Carlo (also called Hamiltonian Monte Carlo, HMC) makes long, directed moves that are still accepted with high probability, by simulating a physical system whose potential energy is \(-\ln p\). It needs the gradient of \(\ln\widetilde p\), and applies to continuous variables.

11.5.1 Dynamical systems

pp. 548–551 Regard \(\mathbf z\) as the position of a particle moving in continuous time \(\tau\), and introduce a momentum (here simply the velocity, unit mass) for each coordinate:

\[ r_i=\frac{\d z_i}{\d\tau}. \tag{11.53} \]

The pair \((\mathbf z,\mathbf r)\) lives in phase space. Write the target as

\[ p(\mathbf z)=\frac1{Z_p}\exp\big(-E(\mathbf z)\big), \tag{11.54} \]

which defines the potential energy \(E(\mathbf z)=-\ln\widetilde p(\mathbf z)\) (any positive density can be written this way). The force is minus the gradient of the potential, and Newton's law says it equals the rate of change of momentum:

\[ \frac{\d r_i}{\d\tau}=-\frac{\partial E(\mathbf z)}{\partial z_i}. \tag{11.55} \]

With the kinetic energy and the Hamiltonian (total energy)

\[ K(\mathbf r)=\frac12\lVert\mathbf r\rVert^2=\frac12\sum_ir_i^2, \tag{11.56} \] \[ H(\mathbf z,\mathbf r)=E(\mathbf z)+K(\mathbf r), \tag{11.57} \]

the dynamics take Hamilton's form

\[ \frac{\d z_i}{\d\tau}=\frac{\partial H}{\partial r_i}, \tag{11.58} \] \[ \frac{\d r_i}{\d\tau}=-\frac{\partial H}{\partial z_i}. \tag{11.59} \]

(Exercise 11.15: \(\partial H/\partial r_i\)\({}=\partial K/\partial r_i\)\({}=r_i\), so (11.58) is (11.53); \(\partial H/\partial z_i=\partial E/\partial z_i\), so (11.59) is (11.55).) Along a trajectory, \(H\) is constant:

\[ \frac{\d H}{\d\tau}=\sum_i\Big\{\frac{\partial H}{\partial z_i}\frac{\d z_i}{\d\tau}+\frac{\partial H}{\partial r_i}\frac{\d r_i}{\d\tau}\Big\}=\sum_i\Big\{\frac{\partial H}{\partial z_i}\frac{\partial H}{\partial r_i}-\frac{\partial H}{\partial r_i}\frac{\partial H}{\partial z_i}\Big\}=0 . \tag{11.60} \]

The second key property, Liouville's theorem, concerns the flow field in phase space,

\[ \mathbf V=\Big(\frac{\d\mathbf z}{\d\tau},\frac{\d\mathbf r}{\d\tau}\Big), \tag{11.61} \]

whose divergence vanishes:

\[ \operatorname{div}\mathbf V=\sum_i\Big\{\frac{\partial}{\partial z_i}\frac{\d z_i}{\d\tau}+\frac{\partial}{\partial r_i}\frac{\d r_i}{\d\tau}\Big\}=\sum_i\Big\{\frac{\partial}{\partial z_i}\frac{\partial H}{\partial r_i}-\frac{\partial}{\partial r_i}\frac{\partial H}{\partial z_i}\Big\}=0 . \tag{11.62} \]

Erratum. In some printings the \(+\) and \(-\) in the second line of (11.62) are swapped; the form above (mixed second derivatives of \(H\) cancelling) is correct.

The joint distribution on phase space that HMC targets is

\[ p(\mathbf z,\mathbf r)=\frac1{Z_H}\exp\big(-H(\mathbf z,\mathbf r)\big). \tag{11.63} \]

Proof

(a) is (11.60).

(b) Write \(\mathbf x=(\mathbf z,\mathbf r)\), \(\dot{\mathbf x}=\mathbf V(\mathbf x)\), and let \(\mathbf J(t)\)\({}=\partial\Phi_t(\mathbf x)/\partial\mathbf x\). Differentiating the equation of motion with respect to the starting point gives \(\dot{\mathbf J}\)\({}=\mathbf D\mathbf V(\Phi_t\mathbf x)\,\mathbf J\), where \(\mathbf D\mathbf V\) is the matrix of partial derivatives of \(\mathbf V\), and \(\mathbf J(0)=\mathbf I\). Jacobi's formula for the derivative of a determinant, \(\frac{\d}{\d t}\det\mathbf J\)\({}=\det\mathbf J\,\Tr(\mathbf J^{-1}\dot{\mathbf J})\), then gives

\[ \frac{\d}{\d t}\det\mathbf J=\det\mathbf J\,\Tr\big(\mathbf J^{-1}\mathbf D\mathbf V\,\mathbf J\big)=\det\mathbf J\,\Tr(\mathbf D\mathbf V)=\det\mathbf J\,\operatorname{div}\mathbf V=0, \]

using the cyclic property of the trace (toolkit) and (11.62). So \(\det\mathbf J(t)=\det\mathbf J(0)=1\) for all \(t\): the image of any region has the same volume as the region.

(c) \(\Phi_t\) is a bijection with inverse \(\Phi_{-t}\) (run the flow backwards). By the change-of-variables formula (11.9) with unit Jacobian, the density of \(\mathbf y=\Phi_t(\mathbf x)\) is \(p(\Phi_{-t}(\mathbf y))\), and by (a) \(H(\Phi_{-t}(\mathbf y))=H(\mathbf y)\), so the density is \(e^{-H(\mathbf y)}/Z_H\), the same as before. This is the book's argument (a small region keeps its volume and its value of \(H\)) made exact.

(d) \(e^{-H}\)\({}=e^{-E(\mathbf z)}e^{-\lVert\mathbf r\rVert^2/2}\) factorises, so \(\mathbf z\) and \(\mathbf r\) are independent with \(p(\mathbf z)\propto e^{-E}\), which is (11.54), and \(p(\mathbf r\given\mathbf z)\)\({}=p(\mathbf r)\)\({}\propto e^{-\lVert\mathbf r\rVert^2/2}\), the standard Gaussian. In particular \(Z_H=Z_p(2\pi)^{D/2}\). ∎

Because \(H\) is constant along trajectories, the dynamics alone are not ergodic: they stay on one energy level. Redrawing \(\mathbf r\) from \(p(\mathbf r\given\mathbf z)\)\({}=\N(\mathbf 0,\mathbf I)\) is a Gibbs step (Section 11.3) that changes \(H\) and leaves (11.63) invariant pp. 550–551. Alternating the two moves and discarding \(\mathbf r\) samples \(p(\mathbf z)\), with the dynamics making large, systematic changes in \(\mathbf z\).

The leapfrog integrator p. 551. The equations must be solved numerically. With step size \(\epsilon\), one leapfrog step is

\[ \widehat r_i(\tau+\epsilon/2)=\widehat r_i(\tau)-\frac{\epsilon}{2}\frac{\partial E}{\partial z_i}\big(\widehat{\mathbf z}(\tau)\big), \tag{11.64} \] \[ \widehat z_i(\tau+\epsilon)=\widehat z_i(\tau)+\epsilon\,\widehat r_i(\tau+\epsilon/2), \tag{11.65} \] \[ \widehat r_i(\tau+\epsilon)=\widehat r_i(\tau+\epsilon/2)-\frac{\epsilon}{2}\frac{\partial E}{\partial z_i}\big(\widehat{\mathbf z}(\tau+\epsilon)\big). \tag{11.66} \]

In a sequence of steps the closing half-step of one and the opening half-step of the next merge into one full momentum step, so positions and momenta "leapfrog" over each other, and only one new gradient is needed per step. Integrating for time \(\tau\) takes \(\tau/\epsilon\) steps.

Proof

(a) Each of (11.64), (11.65), (11.66) changes one block of variables by an amount that depends only on the other block: \((\mathbf z,\mathbf r)\mapsto(\mathbf z,\mathbf r+\mathbf g(\mathbf z))\) or \((\mathbf z,\mathbf r)\mapsto(\mathbf z+\mathbf h(\mathbf r),\mathbf r)\). Such a map is a shear (Figure 11.14 of the book): its Jacobian matrix is block triangular, \(\begin{pmatrix}\mathbf I&\mathbf 0\\\partial\mathbf g/\partial\mathbf z&\mathbf I\end{pmatrix}\) or \(\begin{pmatrix}\mathbf I&\partial\mathbf h/\partial\mathbf r\\\mathbf 0&\mathbf I\end{pmatrix}\), with determinant 1. A composition of maps with unit Jacobian determinant has unit Jacobian determinant (the determinant of a product is the product of the determinants).

(b) Let \((\mathbf z',\mathbf r')\)\({}=\Psi_\epsilon(\mathbf z,\mathbf r)\) with intermediate momentum \(\mathbf r_{1/2}\)\({}=\mathbf r-\frac\epsilon2\nabla E(\mathbf z)\). Apply \(\Psi_{-\epsilon}\) to \((\mathbf z',\mathbf r')\): its first half-step gives \(\mathbf r'+\frac\epsilon2\nabla E(\mathbf z')\)\({}=\mathbf r_{1/2}\) by (11.66); its position step gives \(\mathbf z'-\epsilon\mathbf r_{1/2}=\mathbf z\) by (11.65); its last half-step gives \(\mathbf r_{1/2}\)\({}+\frac\epsilon2\nabla E(\mathbf z)\)\({}=\mathbf r\) by (11.64). So \(\Psi_{-\epsilon}\Psi_\epsilon=\mathrm{id}\), and for \(L\) steps \(\Psi_{-\epsilon}^L\Psi_\epsilon^L=\mathrm{id}\) by undoing the steps one at a time from the last. (Equivalently: flip the momentum, run \(L\) steps forward, flip again.) This is exact in exact arithmetic, whatever \(\epsilon\) is; in floating point it holds to rounding error (in my test with a non-quadratic \(E\), 25 steps forward and back returned to the start within \(5\times10^{-16}\)).

(c) Let \(\mathbf g=\nabla E(\mathbf z)\) and \(\mathbf A=\nabla^2E(\mathbf z)\). The exact flow, by Taylor expansion of (11.53), (11.55): \(\mathbf z(\epsilon)\)\({}=\mathbf z+\epsilon\mathbf r\)\({}-\frac{\epsilon^2}{2}\mathbf g+O(\epsilon^3)\) and \(\mathbf r(\epsilon)\)\({}=\mathbf r-\epsilon\mathbf g\)\({}-\frac{\epsilon^2}{2}\mathbf A\mathbf r\)\({}+O(\epsilon^3)\) (since \(\frac{\d}{\d\tau}\nabla E=\mathbf A\mathbf r\)). Leapfrog: \(\mathbf z'\)\({}=\mathbf z+\epsilon\mathbf r\)\({}-\frac{\epsilon^2}2\mathbf g\) exactly, and \(\mathbf r'\)\({}=\mathbf r-\frac\epsilon2\mathbf g\)\({}-\frac\epsilon2\nabla E(\mathbf z')\)\({}=\mathbf r-\frac\epsilon2\mathbf g\)\({}-\frac\epsilon2\{\mathbf g\)\({}+\epsilon\mathbf A\mathbf r+O(\epsilon^2)\}\)\({}=\mathbf r-\epsilon\mathbf g\)\({}-\frac{\epsilon^2}2\mathbf A\mathbf r\)\({}+O(\epsilon^3)\). The two agree through order \(\epsilon^2\), so the local error is \(O(\epsilon^3)\). Over \(\tau/\epsilon\) steps the local errors accumulate, each propagated by a flow that is Lipschitz over the fixed time interval (Gronwall's inequality), so the global error is \((\tau/\epsilon)\cdot O(\epsilon^3)=O(\epsilon^2)\). Since \(H\) is smooth and exactly conserved by the true flow, \(H\) along the leapfrog path also differs from its initial value by \(O(\epsilon^2)\). In my run on \(E=z^2/2\) up to time 1, halving \(\epsilon\) from 0.1 to 0.05 to 0.025 divided the error in \(z\) by 4.00 and 4.00 (\(3.51\times10^{-4}\), \(8.77\times10^{-5}\), \(2.19\times10^{-5}\)).

(d) Let \(\omega=1/\sigma\). One step is linear, \((z',r')^\T=\mathbf M(z,r)^\T\) with

\[ \mathbf M=\begin{pmatrix}1-\frac{\epsilon^2\omega^2}{2}&\epsilon\\ -\epsilon\omega^2\big(1-\frac{\epsilon^2\omega^2}{4}\big)&1-\frac{\epsilon^2\omega^2}{2}\end{pmatrix}, \]

obtained by composing the three shears. \(\det\mathbf M=1\) (part a) and \(\Tr\mathbf M=2-\epsilon^2\omega^2\). The eigenvalues \(\lambda\) solve \(\lambda^2-(\Tr\mathbf M)\lambda+1=0\). If \(|\Tr\mathbf M|\lt2\), i.e. \(0\lt\epsilon\omega\lt2\), they are complex conjugates of modulus 1 and the iteration is a (sheared) rotation: bounded. If \(\epsilon\omega\gt2\), \(\Tr\mathbf M\lt-2\) and there is a real eigenvalue of modulus \(\gt1\): exponential growth. (In my run, 200 steps with \(\epsilon=1.9\sigma\) stayed bounded, with \(\epsilon=2.1\sigma\) reached \(|z|\approx10^{54}\).) For the conserved quantity, check that \(\mathbf M^\T\mathbf Q\mathbf M=\mathbf Q\) for \(\mathbf Q\)\({}=\diag\big(\omega^2(1-\epsilon^2\omega^2/4),1\big)\), a direct multiplication; so \(\widetilde H=\frac12(z,r)\mathbf Q(z,r)^\T\) is exactly conserved, and \(H-\widetilde H=\frac{\epsilon^2\omega^4}{8}z^2\) is at most \(\dfrac{\epsilon^2\omega^2/4}{1-\epsilon^2\omega^2/4}\,\widetilde H\) (because \(\frac12\omega^2(1-\frac{\epsilon^2\omega^2}4)z^2\)\({}\le\widetilde H\)). In my run with \(\sigma=1\), \(\epsilon=0.3\): \(\widetilde H\) stayed at 0.48875 to 12 digits, while \(H\) oscillated between 0.4888 and 0.5000. ∎

Why gradients help p. 551. The cost of a gradient is usually a small constant multiple of the cost of evaluating \(E\) (reverse-mode differentiation), whatever the dimension \(D\), while it carries \(D\) numbers of information instead of one. Hamiltonian dynamics uses this information to move along the distribution rather than probing at random.

11.5.2 Hybrid Monte Carlo

pp. 552–554 Hybrid Monte Carlo (Duane et al., 1987; Neal, 1996) corrects the leapfrog error with a Metropolis test. From \((\mathbf z,\mathbf r)\), run \(L\) leapfrog steps to a candidate \((\mathbf z^\star,\mathbf r^\star)\) and accept it with probability

\[ \min\big(1,\exp\{H(\mathbf z,\mathbf r)-H(\mathbf z^\star,\mathbf r^\star)\}\big). \tag{11.67} \]

If the integration were exact, \(H\) would not change and every candidate would be accepted. For detailed balance, the integration direction is chosen at random: forwards (\(+\epsilon\)) or backwards (\(-\epsilon\)) with probability \(\frac12\) each.

Proof

The book's version. Both (11.68) and (11.69) equal \(\frac{\delta V}{2Z_H}\min\{e^{-H(\mathcal R)},e^{-H(\mathcal R')}\}\), because \(e^{-a}\min(1,e^{a-b})=\min(e^{-a},e^{-b})\). The volume factor \(\delta V\) is the same in both because the leapfrog map preserves volume (Result 11.25a), and the backward integration from \(\mathcal R'\) lands exactly on \(\mathcal R\) because it inverts the forward one (Result 11.25b). The factor \(\frac12\) is the probability of choosing the right direction. (Svensen and Bishop's official solution splits the cases \(H(\mathcal R)\gt H(\mathcal R')\) and \(\lt\); the identity above does both at once.)

An exact version, which also handles overlapping regions (the book's last remark). Let \(\Psi=\Psi^L_\epsilon\), so \(\Psi^{-1}=\Psi^L_{-\epsilon}\), let \(\pi=e^{-H}/Z_H\), and \(a(\mathbf x,\mathbf y)\)\({}=\min\{1,\pi(\mathbf y)/\pi(\mathbf x)\}\). The accepted-move part of \(\E[g(\mathbf x,\mathbf x')]\) is

\[ \tfrac12\int\pi(\mathbf x)a(\mathbf x,\Psi\mathbf x)\,g(\mathbf x,\Psi\mathbf x)\,\d\mathbf x+\tfrac12\int\pi(\mathbf x)a(\mathbf x,\Psi^{-1}\mathbf x)\,g(\mathbf x,\Psi^{-1}\mathbf x)\,\d\mathbf x . \]

In the second integral substitute \(\mathbf x=\Psi\mathbf y\), with \(\d\mathbf x=\d\mathbf y\) (unit Jacobian): it becomes \(\tfrac12\int\pi(\Psi\mathbf y)a(\Psi\mathbf y,\mathbf y)g(\Psi\mathbf y,\mathbf y)\,\d\mathbf y\). Since \(\pi(\mathbf u)a(\mathbf u,\mathbf v)\)\({}=\min\{\pi(\mathbf u),\pi(\mathbf v)\}\) is symmetric, \(\pi(\Psi\mathbf y)a(\Psi\mathbf y,\mathbf y)\)\({}=\pi(\mathbf y)a(\mathbf y,\Psi\mathbf y)\), so the second integral is the first with the arguments of \(g\) swapped, and by the same substitution the first is the second with the arguments swapped. Hence the accepted part is unchanged when \(g(\mathbf x,\mathbf x')\) is replaced by \(g(\mathbf x',\mathbf x)\). The rejected part has \(\mathbf x'=\mathbf x\) and is trivially symmetric. Taking \(g(\mathbf x,\mathbf x')=h(\mathbf x')\) gives \(\E[h(\mathbf x')]=\E[h(\mathbf x)]\): invariance. The momentum refresh of step 1 is a Gibbs step for (11.63) (Result 11.24d), which also leaves it invariant; marginalising, \(p(\mathbf z)\) is invariant.

Why both directions are needed. With \(+\epsilon\) only, the move \(\mathbf x\to\Psi\mathbf x\) would have no reverse move \(\Psi\mathbf x\to\mathbf x\) (that needs \(\Psi^{-1}\)), and detailed balance would fail. Equivalently one may always integrate forwards and then negate the momentum, which makes the proposal an involution; the negation is harmless because the next step 1 redraws \(\mathbf r\) and \(K(-\mathbf r)=K(\mathbf r)\). ∎

Erratum. In the first printing the arguments of the exponentials inside the minima of (11.68) and (11.69) have the wrong sign; the corrected forms are given above (this is noted in Svensen and Bishop's solution to Exercise 11.17).

Ergodicity p. 553. For some step sizes the leapfrog map returns exactly to its start after \(L\) steps (in the Gaussian case of Result 11.25d this happens when \(L\) times the rotation angle per step is a multiple of \(2\pi\)), and then \(\mathbf z\) never moves however the momentum is refreshed. Drawing \(\epsilon\) at random from a small interval at each iteration rules this out.

11.6 Estimating the partition function

pp. 554–555 Samplers need \(p\) only up to its normalising constant. For model comparison, however, the constant is the quantity of interest. Write

\[ p_E(\mathbf z)=\frac1{Z_E}\exp\big(-E(\mathbf z)\big), \tag{11.71} \]

where \(Z_E\) is the partition function. If \(\mathbf z\) are the parameters \(\bm\theta\) of a model and \(\exp(-E(\bm\theta))\)\({}=p(\mathcal D\given\bm\theta)\,p(\bm\theta)\), then \(p_E\) is the posterior and \(Z_E\)\({}=\int p(\mathcal D\given\bm\theta)p(\bm\theta)\,\d\bm\theta\)\({}=p(\mathcal D)\) is the model evidence (Section 3.4). The ratio of the evidences of two models (the Bayes factor) times the ratio of their prior probabilities is the ratio of their posterior probabilities, so ratios of partition functions are what model comparison needs.

Importance sampling from a distribution \(p_G\) with energy \(G\) gives

\[ \frac{Z_E}{Z_G}=\frac{\sum_{\mathbf z}\exp(-E(\mathbf z))}{\sum_{\mathbf z}\exp(-G(\mathbf z))}=\frac{\sum_{\mathbf z}\exp(-E(\mathbf z)+G(\mathbf z))\exp(-G(\mathbf z))}{\sum_{\mathbf z}\exp(-G(\mathbf z))}=\E_{G}\big[\exp(-E+G)\big]\simeq\frac1L\sum_l\exp\big(-E(\mathbf z^{(l)})+G(\mathbf z^{(l)})\big), \tag{11.72} \]

with \(\mathbf z^{(l)}\) drawn from \(p_G\). If \(Z_G\) is known (a Gaussian, say), this gives \(Z_E\) itself.

Erratum. Some printings omit the factor \(1/L\) on the last line of (11.72), and on the right of (11.73).

A Markov chain as importance distribution p. 555. Analytically tractable \(p_G\) close to \(p_E\) rarely exist. If a Markov chain with transition density \(T(\mathbf z,\mathbf z')\) and invariant distribution \(p_E\) has produced \(\mathbf z^{(1)},\dots,\mathbf z^{(L)}\), one may use

\[ \frac1{Z_G}\exp\big(-G(\mathbf z)\big)=\frac1L\sum_{l=1}^{L}T(\mathbf z^{(l)},\mathbf z) \tag{11.73} \]

in (11.72). This is a mixture of the one-step transition densities from the sampled states. It is normalised (each \(T(\mathbf z^{(l)},\cdot)\) integrates to 1, so \(Z_G=1\) if we take \(\exp(-G)\) equal to the right side), it can be sampled (pick \(l\) at random, then make one transition from \(\mathbf z^{(l)}\)), and, if the \(\mathbf z^{(l)}\) are in equilibrium, its expectation over them is \(\int p_E(\mathbf z')T(\mathbf z',\mathbf z)\,\d\mathbf z'\)\({}=p_E(\mathbf z)\) by invariance (11.39): it is an approximation of \(p_E\) built from the chain. It requires a transition density that can be evaluated, as for a Gibbs sweep with normalised conditionals (a Metropolis transition has a point mass at the current state and is not suitable as it stands).

Chaining p. 555. Ratios can only be estimated between distributions that are close. To reach a complex \(p_M\) from a simple \(p_1\) with known \(Z_1\), insert intermediate distributions \(p_2,\dots,p_{M-1}\) (Neal, 1993; Barber and Bishop, 1997) and write

\[ \frac{Z_M}{Z_1}=\frac{Z_2}{Z_1}\,\frac{Z_3}{Z_2}\cdots\frac{Z_M}{Z_{M-1}}, \tag{11.74} \]

estimating each factor by (11.72) with neighbouring distributions. A convenient path interpolates the energies,

\[ E_\alpha(\mathbf z)=(1-\alpha)E_1(\mathbf z)+\alpha E_M(\mathbf z),\qquad0\le\alpha\le1, \tag{11.75} \]

at a grid of values of \(\alpha\). A single Markov chain can be run through the sequence, moving to the next distribution after some steps, provided it stays close to equilibrium at each stage.

Exercises

All 17 exercises of Chapter 11 (Bishop, 2006, pp. 556–558), restated in my own words, with complete solutions (hidden). Stars give the book's difficulty and "www" marks the exercises with official solutions. Several solutions are compact versions of the results above, which are linked for more detail.

1. Exercise 11.1 (★, www). Show that the estimator \(\widehat f\) of (11.2) has mean \(\E[f]\) and variance (11.3).

   Solution

   Each \(f_l=f(\mathbf z^{(l)})\) has mean \(\E[f]\), so by linearity \(\E[\widehat f]=\frac1L\sum_l\E[f]=\E[f]\). With \(g_l=f_l-\E[f]\), \(\var[\widehat f]=\frac1{L^2}\sum_{l,m}\E[g_lg_m]\); by independence the terms with \(l\ne m\) are \(\E[g_l]\E[g_m]=0\), and the \(L\) diagonal terms are each \(\E[(f-\E f)^2]\). Hence \(\var[\widehat f]=\frac1L\E[(f-\E f)^2]\). Details, and what happens with correlated draws: Result 11.1, Result 11.2.

2. Exercise 11.2 (★). If \(z\) is uniform on \((0,1)\) and \(y=h^{-1}(z)\) with \(h\) given by (11.6), show that \(y\) has density \(p(y)\).

   Solution

   \(h\) is nondecreasing; where \(p\gt0\) it is strictly increasing, so \(\Pr(y\le a)\)\({}=\Pr(h^{-1}(z)\le a)\)\({}=\Pr(z\le h(a))\)\({}=h(a)\), the distribution function of \(p\); differentiating gives \(p(a)\). Equivalently, by (11.5), \(p_y(y)=p_z(h(y))\,|h'(y)|=1\cdot p(y)\). For flat pieces or jumps of \(h\), use the generalised inverse: Result 11.3.

3. Exercise 11.3 (★). Find a transformation \(y=f(z)\) of a uniform \(z\) on \((0,1)\) such that \(y\) has the Cauchy distribution (11.8).

   Solution

   \(h(y)\)\({}=\int_{-\infty}^y\frac{\d t}{\pi(1+t^2)}\)\({}=\frac12+\frac1\pi\arctan y\), so \(f(z)=h^{-1}(z)=\tan\{\pi(z-\frac12)\}\). Check: \(z=\frac12+\frac1\pi\arctan y\) gives \(\d z/\d y=\frac{1}{\pi(1+y^2)}\), which is (11.8) by (11.5). See Filling the gap 11.2.

   ▶ See it in a demo: uniform numbers to Cauchy draws

4. Exercise 11.4 (★★). For \((z_1,z_2)\) uniform on the unit disc and \((y_1,y_2)\) given by (11.10)–(11.11), show that \((y_1,y_2)\) has the density (11.12).

   Solution

   With \(s=r^2\) and \(y_i=z_ig(s)\), \(g(s)=\sqrt{-2\ln s/s}\), the Jacobian matrix is \(g\mathbf I+2g'\mathbf z\mathbf z^\T\), with determinant \(g(g+2sg')=-2/s\). So \(p(y_1,y_2)=\frac1\pi\cdot\frac s2\), and since \(y_1^2+y_2^2=sg^2=-2\ln s\), \(s=e^{-(y_1^2+y_2^2)/2}\) and \(p(y_1,y_2)=\frac1{2\pi}e^{-(y_1^2+y_2^2)/2}\), the product of two standard normal densities. The map is a bijection from the punctured disc onto the punctured plane. Full details: Result 11.4.

   ▶ See it in a demo: Box-Muller

5. Exercise 11.5 (★, www). Let \(\mathbf z\sim\N(\mathbf 0,\mathbf I_D)\) and \(\bm\Sigma=\mathbf L\mathbf L^\T\) (Cholesky). Show that \(\mathbf y=\bm\mu+\mathbf L\mathbf z\) is \(\N(\bm\mu,\bm\Sigma)\).

   Solution

   \(\E[\mathbf y]=\bm\mu\) and \(\cov[\mathbf y]\)\({}=\mathbf L\E[\mathbf z\mathbf z^\T]\mathbf L^\T\)\({}=\mathbf L\mathbf L^\T\)\({}=\bm\Sigma\). A nonsingular affine map of a Gaussian vector is Gaussian: by (11.9) with Jacobian \(|\det\mathbf L|^{-1}=|\bm\Sigma|^{-1/2}\), \(p(\mathbf y)\)\({}=\N(\mathbf L^{-1}(\mathbf y-\bm\mu)\given\mathbf 0,\mathbf I)\,|\bm\Sigma|^{-1/2}\)\({}=\N(\mathbf y\given\bm\mu,\bm\Sigma)\), since \(\mathbf L^{-\T}\mathbf L^{-1}=\bm\Sigma^{-1}\). See Result 11.5.

6. Exercise 11.6 (★★, www). Show that rejection sampling accepts a proposal \(z\) with probability \(\widetilde p(z)/kq(z)\), and that the accepted values have density \(p(z)\).

   Solution

   Given \(z\), \(u\) is uniform on \([0,kq(z)]\), so \(\Pr(u\le\widetilde p(z)\given z)\)\({}=\widetilde p(z)/kq(z)\). The joint density of "propose \(z\), accept" is \(q(z)\widetilde p(z)/kq(z)=\widetilde p(z)/k\); integrating, \(\Pr(\text{accept})=Z_p/k\); dividing, \(p(z\given\text{accept})=\widetilde p(z)/Z_p=p(z)\). (Any valid \(k\) works, not only the smallest.) See Result 11.7.

   ▶ See it in a demo: rejection sampling of a gamma density

7. Exercise 11.7 (★). Show that if \(y\) is uniform, then \(z=b\tan y+c\) has the Cauchy form (11.16).

   Erratum. The printed exercise swaps the roles of \(y\) and \(z\) relative to Section 11.1.2. Also, for the full Cauchy distribution \(y\) must be uniform on \((-\pi/2,\pi/2)\) (or \(y=\pi(u-\frac12)\) with \(u\) uniform on \((0,1)\)); uniform on \([0,1]\) would produce only the values \(c\le z\le c+b\tan1\).

   Solution

   The map \(y\mapsto z=b\tan y+c\) is increasing from \((-\pi/2,\pi/2)\) onto \(\R\), with inverse \(y=\arctan\{(z-c)/b\}\) and \(\d y/\d z=\frac{1/b}{1+(z-c)^2/b^2}\). With \(p(y)=1/\pi\), (11.5) gives

   \[ p(z)=\frac{1}{\pi b}\,\frac{1}{1+(z-c)^2/b^2}, \]

   which is (11.16) with \(k=1/(\pi b)\), the normalised Cauchy density with location \(c\) and scale \(b\). In rejection sampling the constant \(k\) of (11.16) is then rescaled so that \(kq\ge\widetilde p\) (Filling the gap 11.4).

8. Exercise 11.8 (★★). Find the coefficients \(k_i\) of the envelope (11.17) from continuity and normalisation.

   Solution

   On piece \(i\) the unnormalised envelope is \(e^{t_i(z)}\)\({}=\widetilde p(z_i)e^{-\lambda_i(z-z_i)}\), with \(\lambda_i=-\ell'(z_i)\), \(\ell=\ln\widetilde p\). Continuity at \(\widehat z_{i,i+1}\): \(k_i\lambda_ie^{-\lambda_i(\widehat z_{i,i+1}-z_i)}\)\({}=k_{i+1}\lambda_{i+1}e^{-\lambda_{i+1}(\widehat z_{i,i+1}-z_{i+1})}\). Since \(\widehat z_{i,i+1}\) is where the tangents meet, \(\ell(z_i)-\lambda_i(\widehat z-z_i)\)\({}=\ell(z_{i+1})-\lambda_{i+1}(\widehat z-z_{i+1})\), which says exactly that \(k_i\lambda_i/\widetilde p(z_i)\) is the same for all \(i\): \(k_i\lambda_i=\widetilde p(z_i)/C\) for one constant \(C\). Normalisation fixes \(C\):

   \[ C=\sum_i\frac{\widetilde p(z_i)}{\lambda_i}\Big[e^{-\lambda_i(\widehat z_{i-1,i}-z_i)}-e^{-\lambda_i(\widehat z_{i,i+1}-z_i)}\Big],\qquad k_i=\frac{\widetilde p(z_i)}{\lambda_iC}, \]

   with \(\widehat z_{0,1}=-\infty\), \(\widehat z_{n,n+1}=\infty\) (finite because \(\lambda_1\lt0\lt\lambda_n\)). See Result 11.8.

9. Exercise 11.9 (★★). Devise an algorithm for sampling the piecewise exponential distribution (11.17).

   Solution

   (11.17) is a mixture of truncated exponentials. Compute the mass \(M_i\) of each piece (the terms of \(C\) above divided by \(C\)). Draw a piece \(i\) with probability \(M_i\) (inverse distribution function of a discrete distribution). Within the piece \([a,b]\), the distribution function is \(F(z)\)\({}=\{1-e^{-\lambda_i(z-a)}\}/\{1\)\({}-e^{-\lambda_i(b-a)}\}\), so with \(v\) uniform, \(z\)\({}=a-\lambda_i^{-1}\ln[1-v(1-e^{-\lambda_i(b-a)})]\); for the unbounded end pieces, \(z=a-\lambda_n^{-1}\ln v\) (right) and \(z=b+|\lambda_1|^{-1}\ln v\) (left). See Algorithm 11.2 and Filling the gap 11.5.

   ▶ See it in a demo: adaptive rejection sampling

10. Exercise 11.10 (★). For the random walk (11.34)–(11.36), show that \(\E[(z^{(\tau)})^2]=\E[(z^{(\tau-1)})^2]+\frac12\) and hence \(\E[(z^{(\tau)})^2]=\tau/2\).

    Solution

    \(z^{(\tau)}=z^{(\tau-1)}+\xi_\tau\) with \(\xi_\tau\) independent of the past, \(\E\xi_\tau=0\), \(\E\xi_\tau^2=\frac12\). Expanding the square, the cross term has expectation \(2\E[z^{(\tau-1)}]\E[\xi_\tau]=0\), which gives the recursion; with \(z^{(0)}=0\), induction gives \(\tau/2\). See Result 11.15.

    ▶ See it in a demo: distance travelled against \(\sqrt{\text{steps}}\)

11. Exercise 11.11 (★★, www). Show that Gibbs sampling satisfies detailed balance (11.40).

    Solution

    For the update of \(z_k\), \(T(\mathbf z,\mathbf z')\)\({}=p(z_k'\given\mathbf z_{\setminus k})\) when \(\mathbf z'_{\setminus k}=\mathbf z_{\setminus k}\), and 0 otherwise. Then \(p(\mathbf z)T(\mathbf z,\mathbf z')\)\({}=p(\mathbf z_{\setminus k})\,p(z_k\given\mathbf z_{\setminus k})\)\({}\times p(z_k'\given\mathbf z_{\setminus k})\), symmetric in \(z_k\leftrightarrow z_k'\), so it equals \(p(\mathbf z')T(\mathbf z',\mathbf z)\). This holds for each single-variable update and for their random mixture; a fixed-order sweep is a composition, which leaves \(p\) invariant but is reversible only in the forward-backward order. See Result 11.20 and Result 11.18.

12. Exercise 11.12 (★). For the distribution of Figure 11.15 (uniform on two separate regions, one in the lower left and one in the upper right, with no value of \(z_1\) and no value of \(z_2\) in common), is Gibbs sampling ergodic?

    Solution

    No. Suppose the chain is in the lower-left region \(A\). Updating \(z_1\) draws from \(p(z_1\given z_2)\), which is uniform on the horizontal section of the support at height \(z_2\); because the regions share no \(z_2\) values, that section lies entirely in \(A\). Likewise for \(z_2\) updates. So every update stays in \(A\), and the upper-right region \(B\) is never visited: the chain is reducible. It has many invariant distributions (uniform on \(A\), uniform on \(B\), and their mixtures, including the target), and from a start in \(A\) it converges to the uniform distribution on \(A\), not to \(p\). The book's sufficient condition fails because the conditionals vanish on large sets. Remedies: rotate the coordinates (in the coordinates \(z_1+z_2\) and \(z_1-z_2\), a move can cross between the regions if their projections overlap), update \((z_1,z_2)\) as a block, or add a Metropolis-Hastings move that proposes jumps between the regions.

13. Exercise 11.13 (★★). For the graph of Figure 11.16, with \(x\sim\N(x\given\mu,\tau^{-1})\), \(\mu\sim\N(\mu\given\mu_0,s_0)\) and \(\tau\sim\Gam(\tau\given a,b)\) independently, find the conditionals \(p(\mu\given x,\tau)\) and \(p(\tau\given x,\mu)\) needed for Gibbs sampling of \(p(\mu,\tau\given x)\).

    Solution

    Read \(s_0\) as the prior variance of \(\mu\). The joint is \(p(x,\mu,\tau)\)\({}=\N(x\given\mu,\tau^{-1})\N(\mu\given\mu_0,s_0)\Gam(\tau\given a,b)\), and each conditional is proportional to the joint as a function of its own variable (the Markov blanket of \(\mu\) is \(\{x,\tau\}\), and of \(\tau\) is \(\{x,\mu\}\)).

    \(\mu\). Keeping the factors that contain \(\mu\): \(\ln p(\mu\given x,\tau)\)\({}=-\frac\tau2(x-\mu)^2\)\({}-\frac{1}{2s_0}(\mu-\mu_0)^2+\const\)\({}=-\frac12\big(\tau+\frac1{s_0}\big)\mu^2\)\({}+\big(\tau x+\frac{\mu_0}{s_0}\big)\mu+\const\). Completing the square (toolkit):

    \[ p(\mu\given x,\tau)=\N\Big(\mu\;\Big|\;\frac{\tau x+\mu_0/s_0}{\tau+1/s_0},\;\frac{1}{\tau+1/s_0}\Big). \]

    (If \(s_0\) is meant as a precision, replace \(1/s_0\) by \(s_0\).)

    \(\tau\). Keeping the factors that contain \(\tau\): \(\tau^{1/2}e^{-\tau(x-\mu)^2/2}\cdot\tau^{a-1}e^{-b\tau}\)\({}=\tau^{a+1/2-1}e^{-\tau\{b+(x-\mu)^2/2\}}\), the kernel of a gamma density (toolkit):

    \[ p(\tau\given x,\mu)=\Gam\Big(\tau\;\Big|\;a+\frac12,\;b+\frac{(x-\mu)^2}{2}\Big). \]

    Both are standard distributions because the priors are conditionally conjugate; the joint posterior \(p(\mu,\tau\given x)\) is not of standard form (the prior is not the normal-gamma), which is why Gibbs sampling is useful here. With \(N\) observations, \(\tau x\) becomes \(\tau\sum_nx_n\), \(\tau\) becomes \(N\tau\) in the precision, and the gamma parameters become \(a+N/2\) and \(b+\frac12\sum_n(x_n-\mu)^2\).

14. Exercise 11.14 (★). Show that the over-relaxation update (11.50) gives \(z_i'\) with mean \(\mu_i\) and variance \(\sigma_i^2\) when \(z_i\) has mean \(\mu_i\) and variance \(\sigma_i^2\) and \(\nu\) has mean 0 and variance 1.

    Solution

    \(\E[z_i']\)\({}=\mu_i+\alpha(\E z_i-\mu_i)\)\({}+\sigma_i\sqrt{1-\alpha^2}\,\E\nu\)\({}=\mu_i\). With \(z_i\) and \(\nu\) independent, \(\var[z_i']\)\({}=\alpha^2\var[z_i]\)\({}+\sigma_i^2(1-\alpha^2)\var[\nu]\)\({}=\alpha^2\sigma_i^2+\sigma_i^2(1-\alpha^2)\)\({}=\sigma_i^2\). With Gaussian \(z_i\) and \(\nu\), \(z_i'\) is Gaussian too, so its whole conditional distribution is preserved and the update leaves \(p\) invariant (and is reversible): Result 11.22.

    ▶ See it in a demo: ordered over-relaxation

15. Exercise 11.15 (★, www). Using (11.56)–(11.57), show that (11.58) is equivalent to (11.53) and (11.59) to (11.55).

    Solution

    \(H=E(\mathbf z)+\frac12\sum_jr_j^2\), so \(\partial H/\partial r_i=r_i\) (only \(K\) depends on \(\mathbf r\)); thus (11.58), \(\d z_i/\d\tau=\partial H/\partial r_i\), reads \(\d z_i/\d\tau=r_i\), which is (11.53). Similarly \(\partial H/\partial z_i=\partial E/\partial z_i\) (only \(E\) depends on \(\mathbf z\)), so (11.59) reads \(\d r_i/\d\tau=-\partial E/\partial z_i\), which is (11.55).

16. Exercise 11.16 (★). Using (11.56), (11.57) and (11.63), show that \(p(\mathbf r\given\mathbf z)\) is Gaussian.

    Solution

    \(p(\mathbf r\given\mathbf z)\)\({}=p(\mathbf z,\mathbf r)/\int p(\mathbf z,\mathbf r)\,\d\mathbf r\). The factor \(e^{-E(\mathbf z)}/Z_H\) cancels, leaving \(e^{-\lVert\mathbf r\rVert^2/2}/\int e^{-\lVert\mathbf r\rVert^2/2}\d\mathbf r\)\({}=(2\pi)^{-D/2}e^{-\lVert\mathbf r\rVert^2/2}\)\({}=\N(\mathbf r\given\mathbf 0,\mathbf I)\), which does not depend on \(\mathbf z\). See Result 11.24d.

17. Exercise 11.17 (★, www). Verify that (11.68) and (11.69) are equal, so that hybrid Monte Carlo satisfies detailed balance.

    Solution

    For any numbers \(a,b\), \(e^{-a}\min(1,e^{a-b})=\min(e^{-a},e^{-b})\), which is symmetric in \(a\) and \(b\). With \(a=H(\mathcal R)\), \(b=H(\mathcal R')\), (11.68) and (11.69) both equal \(\frac{\delta V}{2Z_H}\min\{e^{-H(\mathcal R)},e^{-H(\mathcal R')}\}\). The equal volumes come from volume preservation of the leapfrog map and the reverse path from its time reversibility, and the factor \(\frac12\) from the random choice of direction. An exact version, valid for overlapping regions: Result 11.26.

    ▶ See it in a demo: energy errors and rejections in hybrid Monte Carlo

Demonstrations for this chapter

Each demonstration runs in the browser, with seeded random numbers so that a class sees the same pictures, and links back to these notes. The ▶ Demo boxes above open them already set to the case under discussion.

The Python code that produces every figure of this page, and checks the chapter's identities numerically, is code/ch11_sampling.py (see the code index).

References

1. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer, New York. Chapter 11, pp. 523–558 (and the official errata).
2. Adler, S. L. (1981). Over-relaxation method for the Monte Carlo evaluation of the partition function for multiquadratic actions. Physical Review D, 23(12), 2901–2904. doi:10.1103/PhysRevD.23.2901.
3. Barber, D. and Bishop, C. M. (1997). Bayesian model comparison by Monte Carlo chaining. In M. Mozer, M. Jordan and T. Petsche (eds), Advances in Neural Information Processing Systems, Vol. 9, 333–339. MIT Press.
4. Box, G. E. P. and Muller, M. E. (1958). A note on the generation of random normal deviates. The Annals of Mathematical Statistics, 29(2), 610–611. doi:10.1214/aoms/1177706645.
5. Celeux, G. and Diebolt, J. (1985). The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem. Computational Statistics Quarterly, 2, 73–82.
6. Devroye, L. (1986). Non-Uniform Random Variate Generation. Springer, New York. doi:10.1007/978-1-4613-8643-8.
7. Duane, S., Kennedy, A. D., Pendleton, B. J. and Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2), 216–222. doi:10.1016/0370-2693(87)91197-X.
8. Fung, R. and Chang, K. C. (1990). Weighting and integrating evidence for stochastic simulation in Bayesian networks. In Uncertainty in Artificial Intelligence, Vol. 5, 208–219. Elsevier.
9. Gelman, A. and Meng, X.-L. (1998). Simulating normalizing constants: from importance sampling to bridge sampling to path sampling. Statistical Science, 13(2), 163–185. doi:10.1214/ss/1028905934.
10. Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6(6), 721–741. doi:10.1109/TPAMI.1984.4767596.
11. Geyer, C. J. (1992). Practical Markov chain Monte Carlo. Statistical Science, 7(4), 473–483. doi:10.1214/ss/1177011137.
12. Gilks, W. R. (1992). Derivative-free adaptive rejection sampling for Gibbs sampling. In J. M. Bernardo, J. O. Berger, A. P. Dawid and A. F. M. Smith (eds), Bayesian Statistics 4. Oxford University Press.
13. Gilks, W. R., Best, N. G. and Tan, K. K. C. (1995). Adaptive rejection Metropolis sampling within Gibbs sampling. Applied Statistics, 44(4), 455–472. doi:10.2307/2986138.
14. Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 41(2), 337–348. doi:10.2307/2347565.
15. Hairer, E., Lubich, C. and Wanner, G. (2006). Geometric Numerical Integration, 2nd edition. Springer, Berlin. doi:10.1007/3-540-30666-8.
16. Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1), 97–109. doi:10.1093/biomet/57.1.97.
17. Henrion, M. (1988). Propagating uncertainty in Bayesian networks by probabilistic logic sampling. In J. F. Lemmer and L. N. Kanal (eds), Uncertainty in Artificial Intelligence, Vol. 2, 149–163. North-Holland. doi:10.1016/B978-0-444-70396-5.50019-4.
18. Jensen, C. S., Kong, A. and Kjærulff, U. (1995). Blocking Gibbs sampling in very large probabilistic expert systems. International Journal of Human-Computer Studies, 42(6), 647–666. doi:10.1006/ijhc.1995.1029.
19. Kong, A. (1992). A note on importance sampling using standardized weights. Technical Report 348, Department of Statistics, University of Chicago.
20. Kong, A., Liu, J. S. and Wong, W. H. (1994). Sequential imputations and Bayesian missing data problems. Journal of the American Statistical Association, 89(425), 278–288. doi:10.1080/01621459.1994.10476469.
21. Levin, D. A., Peres, Y. and Wilmer, E. L. (2017). Markov Chains and Mixing Times, 2nd edition. American Mathematical Society, Providence, RI. doi:10.1090/mbk/107.
22. Liu, J. S. (2001). Monte Carlo Strategies in Scientific Computing. Springer, New York. doi:10.1007/978-0-387-76371-2 (2008 softcover printing).
23. Marsaglia, G. and Bray, T. A. (1964). A convenient method for generating normal variables. SIAM Review, 6(3), 260–264. doi:10.1137/1006063.
24. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H. and Teller, E. (1953). Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6), 1087–1092. doi:10.1063/1.1699114.
25. Metropolis, N. and Ulam, S. (1949). The Monte Carlo method. Journal of the American Statistical Association, 44(247), 335–341. doi:10.1080/01621459.1949.10483310.
26. Neal, R. M. (1993). Probabilistic inference using Markov chain Monte Carlo methods. Technical Report CRG-TR-93-1, Department of Computer Science, University of Toronto.
27. Neal, R. M. (1996). Bayesian Learning for Neural Networks. Lecture Notes in Statistics 118, Springer, New York. doi:10.1007/978-1-4612-0745-0.
28. Neal, R. M. (1999). Suppressing random walks in Markov chain Monte Carlo using ordered overrelaxation. In M. I. Jordan (ed.), Learning in Graphical Models, 205–228. Kluwer, Dordrecht (1998), MIT Press (1999). doi:10.1007/978-94-011-5014-9_8.
29. Neal, R. M. (2001). Annealed importance sampling. Statistics and Computing, 11(2), 125–139. doi:10.1023/A:1008923215028.
30. Neal, R. M. (2003). Slice sampling. The Annals of Statistics, 31(3), 705–767. doi:10.1214/aos/1056562461.
31. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In S. Brooks, A. Gelman, G. L. Jones and X.-L. Meng (eds), Handbook of Markov Chain Monte Carlo, 113–162. Chapman and Hall/CRC. doi:10.1201/b10905-6.
32. Newton, M. A. and Raftery, A. E. (1994). Approximate Bayesian inference with the weighted likelihood bootstrap. Journal of the Royal Statistical Society: Series B, 56(1), 3–48. doi:10.1111/j.2517-6161.1994.tb01956.x.
33. Peskun, P. H. (1973). Optimum Monte-Carlo sampling using Markov chains. Biometrika, 60(3), 607–612. doi:10.1093/biomet/60.3.607.
34. Robert, C. P. and Casella, G. (2004). Monte Carlo Statistical Methods, 2nd edition. Springer, New York. doi:10.1007/978-1-4757-4145-2.
35. Roberts, G. O., Gelman, A. and Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability, 7(1), 110–120. doi:10.1214/aoap/1034625254.
36. Rubin, D. B. (1987). Comment on "The calculation of posterior distributions by data augmentation" by Tanner and Wong: a noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: the SIR algorithm. Journal of the American Statistical Association, 82(398), 543–546. doi:10.2307/2289460.
37. Rubin, D. B. (1988). Using the SIR algorithm to simulate posterior distributions. In J. M. Bernardo, M. H. DeGroot, D. V. Lindley and A. F. M. Smith (eds), Bayesian Statistics 3, 395–402. Oxford University Press.
38. Shachter, R. D. and Peot, M. (1990). Simulation approaches to general probabilistic inference on belief networks. In Uncertainty in Artificial Intelligence, Vol. 5. Elsevier.
39. Tanner, M. A. and Wong, W. H. (1987). The calculation of posterior distributions by data augmentation. Journal of the American Statistical Association, 82(398), 528–540. doi:10.1080/01621459.1987.10478458.
40. Tierney, L. (1994). Markov chains for exploring posterior distributions. The Annals of Statistics, 22(4), 1701–1728. doi:10.1214/aos/1176325750.
41. Wei, G. C. G. and Tanner, M. A. (1990). A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms. Journal of the American Statistical Association, 85(411), 699–704. doi:10.1080/01621459.1990.10474930.

DOIs were checked against Crossref. The numerical statements made in the text (the gamma envelope constants, the importance-sampling table, the random-walk Metropolis and Gibbs autocorrelation times, the leapfrog checks, the chaining variances) come from short Python computations written for these notes; the figures come from the course code.