distribution approaching a multivariate normal in accordance with the
Bernstein–von Mises theorem [177].

#### 12 Informative priors

Finally, the last class of priors are the informative priors the purpose of
which is, unlike Jeffreys and maximum entropy priors which attempt to minimise
the amount of prior information, to include and take into account useful
information for an analysis. They are often highly concentrated in parameter
space and might have been the outcome (i.e. in the form of a posterior
distribution) of a previous experiment or analysis of older data. Their aim is
clearly to inform the analysis and often no attempt is made to restrict the
amount of information provided.

### 3 Making predictions and evaluating models

Tomorrow belongs to those who can hear it coming.
— David Bowie

#### 13 Making predictions

Making predictions is a paramount task for most scientific analyses. Often the
parameters of a model are not observable quantities and we have to rely on
simulated data to assess the validity of our models. In this section, we will
discuss how different kinds of predictive checks can help us avoid various
common pitfalls in Bayesian analyses.

##### 13.1 Prior predictive checks

A very useful practice, that is always recommended, is to check the
predictions of the prior distribution under the specified model [59]. Prior
predictive checks constitute an elegant way of finding out what kind of data
are compatible (i.e. can be described or explained) by our choice of prior and
model. The main benefits of this approach are two. First of all, this can help
diagnose priors that are either too restrictive or too wide. Furthermore,
assuming that the choice of prior distribution is justified, prior predictive
checks can help shield against severe cases of model misspecification in which
no specific set of parameters corresponds to a model that describes the
observed data sufficiently well.

In order to assess whether a particular choice of prior distribution is
appropriate we need a way to produce simulated data that are consistent with
the prior. The Bayesian way of doing this is by sampling the simulated data

$d_{sim}\sim p(d)\,,$ (3.1)

from the prior predictive distribution

$p(d)=\int p(d|\theta)p(\theta)d\theta\,.$ (3.2)

Generating simulated data using the prior predictive distribution in practice
can be done easily by first simulating parameters from the prior distribution

$\theta_{sim}\sim p(\theta)\,,$ (3.3)

and then simulating the data according to the sampling distribution

$d_{sim}\sim p(d|\theta_{sim})\,,$ (3.4)

given the simulated parameters. The simulated pairs $(d_{sim},\theta_{sim})$
constitute samples from the joint distribution

$(d_{sim},\theta_{sim})\sim p(d,\theta)\,,$ (3.5)

and thus

$d_{sim}\sim p(d)\,,$ (3.6)

are simulated from the prior predictive distribution.

##### 13.2 Posterior predictive checks

Similarly to the prior predictive checks, but this time conditioned on the
observed data $d_{obs}$, one can perform posterior predictive checks [59]. The
latter offer a way of measuring whether a model is able to capture aspects of
the data sufficiently well. Just like prior predictive checks that simulate
data consistent with the prior, posterior predictive checks on the other hand
simulate data that are consistent with the posterior distribution.

In practice, the process of generating simulated data

$d_{sim}\sim p(d|d_{obs})\,,$ (3.7)

from the posterior predictive distribution

$p(d|d_{obs})=\int p(d|\theta)p(\theta|d_{obs})d\theta\,,$ (3.8)

starts by simulating parameters from the posterior distribution

$\theta_{sim}\sim p(\theta|d_{obs})\,.$ (3.9)

It is important to remind the reader that this last step, unless conjugate
priors are used, is highly non–trivial and often requires advanced
computational algorithms that are the subject of the next chapter. For now, it
suffices to understand that the simulation of parameters as described by the
relation 3.9 is possible although generally difficult, requiring careful
steps. The last step is to generate the simulated data according to sampling
distribution

$d_{sim}\sim p(d|\theta_{sim})\,,$ (3.10)

given the simulated parameters. $d_{sim}$ then constitute samples from the
posterior predictive distribution.

#### 14 Evaluating and comparing models

A key goal of science is to determine which model under consideration better
accounts for the observed data. As we will discover shortly, this is generally
done by assessing the predictive power of different models. From a Bayesian
perspective, there are two ways one can approach this subject. The first uses
Bayes factors and compares models based on their prior predictive performance,
meaning their capacity to explain the observed data using only the information
encoded in the prior distribution. On the other hand, the second approach uses
the notion of cross–validation in order to compare models based on their
posterior predictive performance, meaning their ability to make predictions
for out–of–sample data, meaning, future or unseen data, using what we learned
from the observed data.

In the prior predictive approach, the main quantity that goes into the
calculation of the Bayes factor is the prior predictive probability
$p(d|\mathcal{M}_{i})$ of the observed data $d$ given a model
$\mathcal{M}_{i}$, also known as the marginal likelihood or the model evidence
[92, 74]. Naturally, the prior predictive approach is sensitive to the choice
of priors. On the other hand, in the posterior predictive approach, we compute
the posterior predictive probability of some subset of the observed data given
the rest of the data. Typically, cross–validation means that this process is
repeated several times, trying to predict different subsets of data, until the
entire data set is assessed as held–out data.

##### 14.1 Bayes factors

The probability of a model $\mathcal{M}_{i}$ given the data $d$ can be
computed using Bayes’ theorem

$p(\mathcal{M}_{i}|d)=\frac{p(d|\mathcal{M}_{i})p(\mathcal{M}_{i})}{p(d)}\,,$
(3.11)

where $p(d|\mathcal{M}_{i})$ is the probability of the data given the model,
$p(\mathcal{M}_{i})$ is the prior probability of the model, and $p(d)$ is the
prior predictive probability of the data. We can compare two models,
$\mathcal{M}_{i}$ and $\mathcal{M}_{j}$, by computing their odds ratio

$\frac{p(\mathcal{M}_{i}|d)}{p(\mathcal{M}_{j}|d)}=\frac{p(d|\mathcal{M}_{i})p(\mathcal{M}_{i})}{p(d|\mathcal{M}_{j})p(\mathcal{M}_{j})}$
(3.12)

The ratio $BF_{ij}=p(d|\mathcal{M}_{i})/p(d|\mathcal{M}_{j})$ in the above
expression is called the Bayes factor. Once the model priors
$p(\mathcal{M}_{i})$ and $p(\mathcal{M}_{j})$ are specified, model comparison
using Bayes factors amounts to the calculation of the model evidences
$p(d|\mathcal{M}_{i})$ and $p(d|\mathcal{M}_{j})$.

Figure 3.1: Prior predictive distributions $p(d|\mathcal{M}_{i})$ and
$p(d|\mathcal{M}_{j})$ for models ${M}_{i}$ and ${M}_{j}$ respectively. The
dashed line, that intersects both distributions, corresponds to the actual
observed data. The Bayes factor is simply the ration between the values at the
two points of intersection, in this case favouring ${M}_{j}$ over ${M}_{i}$.
It is clear that for other realisations of the actual observed data (e.g. on
the right part of the data vector) the other model would be favoured.

Despite the apparent simplicity of model comparison using Bayes factors,
caution must be exercised when applying the method in real analyses. There are
three main reasons for this warning, all of which are sometimes overlooked in
practice leading to catastrophic results.

The first reason has to do with the computational difficulty of estimating the
model evidence $p(d|\mathcal{M})$, particularly in problems with many
parameters. In fact, as we will see in detail in the next part of this thesis,
a large collection of methods have been designed with the sole purpose of
estimating the model evidence. Therefore, the practitioner has to be familiar
with the range of applicability of each method as well as their intrinsic
limitations when deciding which technique to use.

The second reason, equally important with the first, is the sensitivity of the
model evidence to the choice of prior distribution. This sensitivity is
apparent if we just notice that the model evidence is simply the prior
predictive distribution,

$p(d|\mathcal{M})=\int p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})d\theta\,,$
(3.13)

evaluated at the observed data $d$. However, we argue that this sensitivity is
not a weakness of the method as it is often portrayed, but a strength that
needs to be properly understood.

Figure 3.2: The characteristic width $\delta\theta$ of the likelihood function
$p(d|\theta,\mathcal{M}_{1})$ and $\Delta\theta$ of the prior distribution.

In order to understand the effects of the choice of priors on the Bayes
factor, let us consider a simple example. Imagine that we have to compare two
models, $\mathcal{M}_{1}$ with a single scalar parameter $\theta$ and
$\mathcal{M}_{0}$ with no free parameters. Furthermore, let us assume that
$\mathcal{M}_{0}$ is nested in $\mathcal{M}_{1}$, meaning that the more
complex model, $\mathcal{M}_{1}$, reduces to the simpler one,
$\mathcal{M}_{0}$, for a specific parameter value, $\theta=\theta_{0}$.

Let us now assume that the prior on parameter $\theta$ is flat or uniform,
such that,

$p(\theta|\mathcal{M}_{1})=\frac{1}{\Delta\theta}\,,$ (3.14)

and that the likelihood function is sharply peaked around a value $\theta_{1}$
such that,

$\int
p(d|\theta,\mathcal{M}_{1})d\theta=p(d|\theta_{1},\mathcal{M}_{1})\times\delta\theta\,,$
(3.15)

where $\delta\theta$ is its characteristic width. It is easy to show that for
the case of Gaussian likelihood, centred around $\theta_{1}$, the
characteristic width is simply $\delta\theta=\sqrt{2\pi}\sigma$, where
$\sigma$ is the standard deviation.

The model evidence of $\mathcal{M}_{1}$ is then simply,

$\begin{split}p(d|\mathcal{M}_{1})&=\int
p(d|\theta,\mathcal{M}_{1})p(\theta|\mathcal{M}_{1})d\theta\\\
&=p(d|\theta_{1},\mathcal{M}_{1})\times\frac{\delta\theta}{\Delta\theta}\,.\end{split}$
(3.16)

Since the model $\mathcal{M}_{0}$ has no free parameters, no integration is
required and its model evidence $p(d|\mathcal{M}_{0})$ is simply the
likelihood function of $\mathcal{M}_{1}$ evaluated at $\theta=\theta_{0}$, or,

$p(d|\mathcal{M}_{0})=p(d|\theta_{0},\mathcal{M}_{1})\,.$ (3.17)

Therefore, the Bayes factor is,

$B_{10}=\frac{p(d|\mathcal{M}_{1})}{p(d|\mathcal{M}_{0})}=\frac{p(d|\theta_{1},\mathcal{M}_{1})}{p(d|\theta_{0},\mathcal{M}_{1})}\times\frac{\delta\theta}{\Delta\theta}\,.$
(3.18)

The first term in equation 3.18 is the likelihood ratio that always favours
the most complex model $\mathcal{M}_{1}$ since it contains $\mathcal{M}_{0}$
as a special case. In other words, the first term is always greater than one
as the most complex model can always fit the data better than the simpler one.

On the other hand, the second term in equation 3.18 that consists of the ratio
of the likelihood width $\delta\theta$ to the prior width $\Delta\theta$
penalises the most complex model $\mathcal{M}_{1}$, since
$\delta\theta<\Delta\theta$, for any “wasted” regions of parameter space that
are ruled out by the data. This term quantifies the so–called principle of
parsimony or Occam’s razor as it most commonly known. The principle, often
attributed to William of Ockham, states that “entities should not be
multiplied beyond necessity”, meaning that between competing models or
hypotheses the simplest one is often preferred. Therefore, the Bayes factor
will only favour the most complex model (i.e. $\mathcal{M}_{1}$) only if the
likelihood ratio is large enough to overcome the penalty introduced by Occam’s
razor. This intrinsic property of Bayes factors to prefer simpler models, that
originates directly from the reliance to the prior distributions, is what
makes them so useful in practice.

Now that we understand how sensitive the Bayes factor is to the choice of
priors we can discuss some ways that we can shield our analyses against
potential problems. First and foremost, Bayesian model comparison can be
performed only when proper priors are used. By that we mean that improper
priors such as uniform/flat priors ranging from $-\infty$ to $+\infty$ are not
acceptable. Only prior distributions that can be integrated and normalised to
unity are proper in this sense. However, the use of proper priors is not
enough, the choice of priors needs to be well–justified too. Priors that are
not defined using a principled process (e.g. MaxEnt, Jeffreys, etc.), and
sometimes even those that do, can lead to significant deviations in the value
of a Bayes factor. For this reason, we urge caution not to over–emphasise the
significance of, and instead mostly neglect Bayes factors of $\mathcal{O}(1)$.

The third, and final in our list of reasons, has to do with the open–ended
nature of the task model comparison. In particular, model comparison often
takes place in the context of a finite set of possible models under
investigation with no guarantee whatsoever that one of those models captures
perfectly, or even sufficiently, the true data generating process. In that
sense, in almost all cases, inference takes place under conditions of model
misspecification. This brings to mind the saying by Box, that “all models are
wrong, but some are useful”. The fact that the value of a Bayes factor might
seem to favour one model over another does not mean that the first model is
“correct”, only that it is better than the second. Both models might be far
from the true data generating process and the Bayes factor will offer
generally no indication of that.

##### 14.2 Cross–validation

For a model to be useful in practice it must be able to make accurate
predictions regarding unseen data. The generalisation uncertainty of a model
is often quantified using some measure of the out–of–sample predictive
accuracy. A commonly used scoring rule for the out–of–sample predictive
accuracy for $n$ data points is the expected log–pointwise predictive density,

$\mathrm{ELPD}=\sum_{i=1}^{n}\int p_{t}(d_{i})\log p(d_{i}|d_{obs})dd_{i}\,,$
(3.19)

where $p_{t}(d_{i})$ is the probability density of the true data generative
process which is in general unknown and $p(d_{i}|d_{obs})$ is the posterior
predictive density.

Another useful quantity is the log–pointwise predictive density,

$\mathrm{LPD}=\sum_{i=1}^{n}\log p(d_{i}|d_{obs})=\sum_{i=1}^{n}\log\int
p(d_{i}|\theta)p(\theta|d_{obs})d\theta\,.$ (3.20)

LPD of the observed data $d_{obs}$ is an overestimate of ELPD. We can compute
LPD in practice as,

$\hat{\mathrm{LPD}}=\sum_{i=1}^{n}\log\left(\frac{1}{J}\sum_{j=1}^{J}p(d_{i}|\theta_{j})\right)\,,$
(3.21)

where $\theta_{j}\sim p(\theta|d_{obs})$ are samples from the posterior
distribution.

###### 14.2.1 Leave–one–out cross–validation

The term cross–validation refers to the practice of estimating the
out–of–sample predictive accuracy of a model. In general, the method requires
running the analysis multiple times, each time excluding a different portion
of the data. The excluded part of the data is then used in order to assess the
predictive accuracy of the model. Once the whole dataset is covered, the total
accuracy is computed as the average accuracy over all runs,

$\mathrm{ELPD}_{\mathrm{LOO}}=\sum_{i=1}^{n}\log p(d_{i}|d_{-i})\,,$ (3.22)

where,

$p(d_{i}|d_{-i})=\int p(d_{i}|\theta)p(\theta|d_{-i})d\theta\,,$ (3.23)

is the leave–one–out predictive density given the data without the $i$–th
datapoint [55, 12, 70].

Assuming that the $n$ datapoints are conditionally independent in the data
generative model, then we can approximate equation 3.23 using draws from the
posterior $\theta_{j}\sim p(\theta|d_{obs})$ and importance weights [56],

$w_{ij}=\frac{1}{p(d_{i}|\theta_{j})}\propto\frac{p(\theta_{j}|d_{-i})}{p(\theta_{j}|d_{obs})}\,,$
(3.24)

leading to the importance sampling leave–one–out predictive density,

$p(d_{i}|d_{-i})=\frac{\sum_{j=1}^{J}w_{ij}p(d_{i}|\theta_{j})}{\sum_{j=1}^{J}w_{ij}}=\frac{1}{\frac{1}{J}\sum_{j=1}^{J}[p(d_{i}|\theta_{j})]^{-1}}\,.$
(3.25)

However the posterior $p(\theta|d_{obs})$ is likely to have a smaller variance
than then $p(\theta|d_{-i})$ distributions leading to insufficient overlap
between their typical sets and high–variance importance weights. [90] showed
that truncating the importance weights,

$\tilde{w}_{ij}=\min\left(w_{ij},\sqrt{J}\,\bar{w}_{i}\right)\,,$ (3.26)

where

$\bar{w}_{i}=\frac{1}{J}\sum_{j=1}^{J}w_{ij}\,,$ (3.27)

leads to provable finite–variance weights at the cost of introducing bias.
[179] proposed instead to fit a generalised Pareto distribution to the upper
tail of the importance weights, in order to smooth the weights, leading to
improved estimates.

###### 14.2.2 WAIC

The Watanabe–Akaike or widely applicable information criterion (WAIC) [185]
offers a different way to approximate ELPD and is defined as,

$\hat{\mathrm{ELPD}}_{\mathrm{WAIC}}=\hat{\mathrm{LPD}}-\hat{p}_{\mathrm{WAIC}}\,,$
(3.28)

where,

$\hat{p}_{\mathrm{WAIC}}=\sum_{i=1}^{n}\mathrm{Var}_{\theta\sim
p(\theta|d_{obs})}[\log p(d_{i}|\theta)]\,,$ (3.29)

is the estimated effective number of parameters expressed as the posterior
variance of the log predictive density of each datapoint. Equation 3.29 can be
computed using posterior samples.

##### 14.3 Model averaging

Standard practice ignores model uncertainty and instead focuses on the most
probable models as deduced by their Bayes factors. This approach leads to
over–confident estimates and ignores the fact that often more than one model
can describe the data sufficiently. There is, however, a different approach
that we can follow in order to deal with the model uncertainty and properly
account for the plethora of plausible models, called Bayesian model averaging
[119].

Let $\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{M}$ be a set of $M$
models with posterior model probabilities
$p(\mathcal{M}_{1}|d),\,p(\mathcal{M}_{2}|d),\,\dots,\,p(\mathcal{M}_{M}|d)$
and posterior distributions $p(\theta|d,\mathcal{M}_{1}),\,\allowbreak
p(\theta|d,\mathcal{M}_{2}),\,\allowbreak\dots,\,\allowbreak
p(\theta|d,\mathcal{M}_{M})$ respectively. Then Bayesian model averaging
relies on the marginal posterior density,

$p(\theta|d)=\sum_{i=1}^{M}p(\theta|d,\mathcal{M}_{i})p(\mathcal{M}_{i}|d)\,,$
(3.30)

which is no longer conditioned on a model.

Moreover, predictions can be made by averaging over all models, weighted
proportional to their posterior model probabilities, thereby incorporating
model uncertainty using the marginal posterior predictive density,

$p(d|d_{\mathrm{obs}})=\sum_{i=1}^{M}p(d|d_{\mathrm{obs}},\mathcal{M}_{i})p(\mathcal{M}_{i}|d_{\mathrm{obs}})\,,$
(3.31)

where $d_{\mathrm{obs}}$ are the available observed data and $d$ are the new
predicted data. [118] note that averaging over all models in this fashion
leads to higher predictive accuracy than using any single model individually.

##  Part II Bayesian Computation

### 4 Principles of Bayesian Computation

Anyone who considers arithmetical methods of producing
random digits is, of course, in a state of sin.
— John von Neumann

This chapter introduces the various methods that are used in practice in order
to tackle the computational challenges of Bayesian analyses. We begin this
journey by discussing some fundamental ideas about the geometry of
high–dimensional probability distributions, while gradually introducing the
concepts and algorithms that constitute the modern mathematical machinery of
Bayesian computation.

#### 15 Expectation values

Probability theory teaches us the only well defined way to extract information
from probability distributions is through expectation values. By this term, we
mean high–dimensional integrals of the form

$\mathbb{E}_{p}[f]=\int f(\theta)p(\theta)d\theta\,,$ (4.1)

where $p(\theta)$ is the probability density function that often corresponds
to the posterior density for problems of scientific inference, $\theta$
signifies the parameters of the distribution, and $f(\theta)$ is the function
that we aim to integrate. In this sense, an expectation value of a function
$f(\theta)$ over a probability distribution $p(\theta$) is technically a
functional of the product of the function and the probability density.

To see why expectation values hold such a central role in scientific parameter
inference, let us discuss a few characteristic and common examples that a
scientist often has to compute.

  * $\bullet$

Mean value – Perhaps the most commonly computed expectation value is the mean
value. This can be calculated by choosing the function to be
$f(\theta)=\theta$, the expectation value then reduces to

$\mu\equiv\mathbb{E}_{p}[\theta]=\int\theta\,p(\theta)d\theta\,.$ (4.2)

  * $\bullet$

Variance – One might also want to compute higher moments of the probability
distribution. The first moment is the variance that corresponds to the
following expectation value

$\sigma^{2}\equiv\mathbb{E}_{p}\left[(\theta-\mu)^{2}\right]=\int(\theta-\mu)^{2}\,p(\theta)d\theta\,.$
(4.3)

  * $\bullet$

Marginal distributions – Even marginal distribution can be thought of as
expectation values. In this case, the function $f$ corresponds to a
conditional distribution, for instance

$p(\phi)\equiv\mathbb{E}_{p}\left[p(\phi|\theta)\right]=\int
p(\phi|\theta)p(\theta)d\theta\,.$ (4.4)

#### 16 Quadrature and uniform grids

Figure 4.1: Uniform grid approximation of 1–dimensional probability
distribution.

Having discussed a number of examples of commonly used expectation values, we
can now turn our attention to the methods that are used for their computation.
As we mentioned before, these expectation values are defined as
high–dimensional integrals. As those integrals are not generally tractable
analytically, one might attempt to approximate their value by means of a
Riemann sum over a discreet grid of $n$ points:

$\mathbb{E}_{p}[f]=\int
f(\theta)p(\theta)d\theta\approx\sum_{i=1}^{n}f(\theta_{i})p(\theta_{i})\Delta\theta_{i}\,,$
(4.5)

where

$\Delta\theta_{i}=\theta_{j+1}-\theta_{j}\,,$ (4.6)

is simply the interval between two subsequent points, $\theta_{j}$ and
$\theta_{j+1}$ on the underlying grid, and

$\theta_{i}=\frac{\theta_{j+1}+\theta_{j}}{2}\,,$ (4.7)

is just the mid–point between $\theta_{j}$ and $\theta_{j+1}$.

In principle, this idea can be extended to high dimensions by replacing the
1–dimensional intervals $\Delta\theta_{i}$ with D–dimensional hypercubes.
Figure 4.2 shows one such example for a 2–dimensional probability
distribution. However, as the number of dimensions increases, one immediately
has to face a significant difficulty, the _curse of dimensionality_ [8].
Already in 2 dimensions we require $n^{2}$ grid points to approximate the
distribution. As it turns out, the number of grid points required for the
evaluation of the Riemann sum increases exponentially with the number of
dimensions, rendering this method of computing expectation values unusable for
$D>3$. Overcoming the difficulties imposed by the _curse of dimensionality_ is
one of the key goals of probabilistic computing.

Figure 4.2: Uniform grid approximation of 2–dimensional probability
distribution.

In order to reduce the computational cost of estimating expectation values in
high dimensions, we need to find a way to focus our effort and computation
only on those regions of parameter space that are relevant for the integral
that we aim to evaluate. One simple idea would be to remove any points of the
grid that the integrand $f(\theta)p(\theta)$ is very close to zero. Applying
this technique would certainly reduce the total computational cost since only
a few grid–cells have a non–negligible value of $f(\theta)p(\theta)$ as shown
in Figure 4.3. The problem that we face however is that by focusing our
attention on $f(\theta)p(\theta)$ we ignore a key factor in the estimation of
any expectation value, the volume.

Figure 4.3: Uniform grid approximation of 2–dimensional probability
distribution with highlighted the grid–cells that actually contribute to the
calculation of an expectation value.

#### 17 Geometry of high–dimensional spaces

The concepts of volume and distance in high–dimensional spaces defy our
everyday intuition in ways that matter for the computation of expectation
values. To understand this, we will go through an example that illustrates
these peculiar effects.

Let us assume that we inscribe a circle of radius $R$ inside a square of side
$2R$. We are interested in computing the area of the circle as a fraction of
that of the square. We can get to the result easily using basic geometry, in
particular, the ratio of the two areas is

$\frac{A_{\textrm{circle}}}{A_{\textrm{square}}}=\frac{\pi
R^{2}}{(2R)^{2}}=\frac{\pi}{4}\,.$ (4.8)

We can now extend the same problem into three dimensions, in which we have a
sphere of radius $R$ inscribed in a cube of side $2R$. The ratio of the volume
of the sphere to the volume of the cube is simply

$\frac{V_{\textrm{sphere}}}{V_{\textrm{cube}}}=\frac{\frac{4}{3}\pi
R^{3}}{(2R)^{3}}=\frac{\pi}{6}\,.$ (4.9)

By comparing equations 4.8 and 4.9 one realises that the volume ratio has
decreased going from $2$ dimensions to $3$. We will now show that this result
in fact holds for any number of dimensions $D$. In $D$ dimensions, the volume
of a hyper–sphere is given by

$V_{\textrm{sphere}}=\frac{\pi^{D/2}}{\Gamma\left(\frac{D}{2}+1\right)}R^{D}\,,$
(4.10)

where $\Gamma$ is Euler’s gamma function which extends the factorial operation
to non–integer arguments and satisfies

$\Gamma(D)=(D-1)!\,,$ (4.11)

for positive integer $D$, and

$\Gamma\left(D+\frac{1}{2}\right)=\left(D-\frac{1}{2}\right)\times\left(D-\frac{3}{2}\right)\times\dots\times\frac{1}{2}\times\pi^{1/2}\,,$
(4.12)

for non–negative integer $D$.

The volume of a hyper–cube in $D$ dimensions is simply

$V_{\textrm{cube}}=(2R)^{D}\,.$ (4.13)

Taking the ratio of the terms of equations 4.10 and 4.13 yields

$\frac{V_{\textrm{sphere}}}{V_{\textrm{cube}}}=\frac{\pi^{D/2}}{2^{D}\Gamma\left(\frac{D}{2}+1\right)}\,.$
(4.14)

Figure 4.4 shows the ratio of the volume of a hypersphere to a hypercube as a
function of the number of dimensions $D$. As the number of dimensions $D$
increases, the volume ratio of equation 4.14 asymptotically approaches $0$.
This means that in high dimensions, almost all of the volume of a hypercube is
concentrated in the corners.

Figure 4.4: The ratio of the volume of a hyper–sphere of radius $R$ to the
volume of a hyper–cube of edge size $2R$ as function of the number of
dimensions $D$.

#### 18 Concentration of measure

As we will see shortly, the strange behaviour of volume is of paramount
importance in the calculation of expectation values over probability
distributions. To understand this one need to think not about the probability
density but instead about the probability mass. When evaluating an expectation
value, not all regions of parameter space are contributing equally to the
value of the integral. In fact, the contribution from some regions of
parameter space dominates the calculation. We only need to take a look into
the form of the expectation value integral to notice that is essentially the
product of three terms that contributes. These terms are the function
$f(\theta)$, the probability density function $p(\theta)$, as well as the
differential volume element $d\theta$. In other words, it is the product of
these three terms and their dependence on $\theta$ that determines the value
of the integral. Assuming that the function $f(\theta)$ is well behaved, we
can ignore its presence for a while.

For the sake of simplicity let us assume that the probability distribution is
characterised by an D–dimensional Gaussian probability density function

$p(\theta)=\det(2\pi\Sigma)^{-\frac{1}{2}}\exp{\left[-\frac{1}{2}(\theta-\mu)^{T}\Sigma^{-1}(\theta-\mu)\right]}\,,$
(4.15)

where $\mu$ is the mean and $\Sigma$ is the covariance matrix of the
distribution. Without loss of generality let us also assume that the density
is centred at zero (i.e. $\mu=0$) and the covariance matrix diagonal with the
elements of its diagonal equal to $\sigma^{2}$, meaning that equation 4.15
simplifies into

$p(\theta)=\frac{1}{\sqrt{(2\pi)^{D}}\sigma^{D}}\exp{\left(-\frac{|\theta|^{2}}{2\sigma^{2}}\right)}\,.$
(4.16)

Assuming further spherical coordinates, the density only depends on the
magnitude $r$ of the $\theta$ parameter vector

$p(r)=\frac{1}{(2\pi)^{\frac{D}{2}}\sigma^{D}}\exp{\left(-\frac{r^{2}}{2\sigma^{2}}\right)}\,.$
(4.17)

Keep in mind that $p(r)$ is not a probability density function of the
magnitude $r=|\theta|$, but a D–dimensional density of $\theta$.

Let us now turn our attention to the differential volume element $dV$.
Differentiating equation 4.10 that provides the volume of the hyper–sphere we
get

$dV=\frac{D\pi^{D/2}}{\Gamma\left(\frac{D}{2}+1\right)}r^{D-1}dr\,.$ (4.18)

Figure 4.5: Scaling of differential volume with the number of dimensions as a
function of distance.

The differential probability mass $dm(r)$ is then just the product of $p(r)$
and $dV$ given by equations 4.17 and 4.18 respectively

$dm(r)=\frac{Dr^{D-1}}{\Gamma\left(\frac{D}{2}+1\right)2^{\frac{D}{2}}\sigma^{D}}\exp{\left(-\frac{r^{2}}{2\sigma^{2}}\right)}dr\,.$
(4.19)

The differential mass $dm(r)$ has a clear physical meaning, that of the
probability mass enclosed in a hyper–spherical shell of radius $r$ and width
$dr$. The probability mass differential $dm(r)$ peaks (i.e. is maximised) at
the typical radius

$r_{\textrm{peak}}=\sqrt{D-1}\sigma\,.$ (4.20)

Equation 4.20 indicates that while in 1–D the probability mass peaks at
$r_{\textrm{peak}}=0$, in higher dimensions this is not the case. As the
number of dimensions increases the radius in which the probability mass peaks
moves to greater distances. Table 1 shows the typical radius of the
probability mass for a different number of dimensions for our problem. This is
a direct consequence of the rapid increase of the differential volume for
large $r$ values.

Figure 4.6: Scaling of differential probability mass with the number of dimensions as a function of distance. Table 1: The typical radius $r_{\textrm{peak}}$ as function of the number of dimensions $D$. [0.75pt] Number of dimensions $D$ | Typical radius $r_{\textrm{peak}}$
---|---
[0.5pt] $1$ | $0$
$2$ | $1\sigma$
$5$ | $2\sigma$
$10$ | $3\sigma$
$17$ | $4\sigma$
$26$ | $5\sigma$
[0.75pt] |

In general, we expect the probability mass to form a hyper–shell of mean
radius

$r_{\textrm{mean}}\equiv\mathbb{E}_{p}[r]=\int_{0}^{+\infty}rdm(r)$ (4.21)

and width (i.e. standard deviation)

$\Delta
r\equiv\sqrt{\mathbb{E}_{p}[(r-r_{\textrm{mean}})^{2}]}=\sqrt{\int_{0}^{+\infty}(r-r_{\textrm{mean}})^{2}dm(r)}$
(4.22)

Figure 4.7: Scaling of differential probability mass with the number of
dimensions as a function of distance normalised by the square root of the
number of dimensions.

#### 19 Typical set

Figure 4.8: Illustration of the typical set as the region in parameter space
that the product of probability density and differential volume is
non–negligible.

The qualitative conclusions of the previous section are general and hold for
any continuous probability distribution. The probability mass does not
concentrate close to the mode where the probability density is high as there
is not sufficient volume there. On the other hand, it does not concentrate on
large distances because the density vanishes. Instead, it compromises on some
region of intermediate distance surrounding the mode, as shown in Figure 4.8.
This region is called the typical set, and has the form of a high–dimensional
thin hyper–shell surrounding the mode as shown in Figure 4.9. In high
dimensions, the typical set exhibits the effect of concentration of measure
[109] illustrated in Figure 4.7.

Figure 4.9: Illustration of the typical set as a thin hyper–shell surrounding
the mode of the probability distribution.

The concept of the typical set is not only important for properly
understanding probability distributions, but also for developing new
computational methods. The notion of the typical set is originally borrowed
from the field of information theory, in which one of the main tasks is to
compress and encode a message with as few words as possible. In probability
theory, the typical set defines the most efficient way to compress a
probability distribution by focusing on a limited region of parameter space.
As we will see in the next sections, the task of developing powerful and
effective computational methods comes down to how efficiently we can locate
and approximate the typical set of a probability distribution.

#### 20 Laplace approximation

Before we move on to stochastic estimators of expectation values let us first
discuss another simple deterministic method, called Laplace approximation,
that, unlike quadrature in a uniform grid, can extend to higher dimensions
[173]. The Laplace approximation makes a very strong assumption about the
target probability distribution. In particular, it assumes that it can be
sufficiently described by a Gaussian probability density, similar to equation
4.15. The mean of the Gaussian density is determined at the point of the mode
of the target density

$\mu=\underset{\theta}{\mathrm{arg\,max}}\leavevmode\nobreak\ p(\theta)\,,$
(4.23)

and the precision matrix $\Sigma^{-1}$ (i.e. inverse of the covariance matrix
$\Sigma$) is given by the second-order derivatives of the negative logarithm
of the target probability density function evaluated at the mode,

$\left(\Sigma^{-1}\right)_{ij}=-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\log
p(\theta)\bigg{|}_{\theta=\mu}\,.$ (4.24)

Figure 4.10: Illustration of the Laplace approximation to a skewed probability
density.

The reasoning behind this approach is quite simple, one effectively performs a
Taylor expansion of the logarithm of the density, up to second order, around
the maximum a posteriori point $\mu$,

$\log p_{L}(\theta)=\log
p(\theta=\mu)-\frac{1}{2}(\theta-\mu)^{T}\Sigma^{-1}(\theta-\mu)+\dots\,,$
(4.25)

where the first order term simply vanishes because we evaluate the expansion
around the maximum. For this reason, this very common method is often called
the saddle–point approximation. Expectation values can then be determined
using the Gaussian density $p_{L}(\theta)=\mathcal{N}(\theta|\mu,\Sigma)$ in
place of the target density $p(\theta)$ in the formula for the expectation
value 4.1,

$\mathbb{E}_{p_{L}}[f]=\int f(\theta)\mathcal{N}(\theta|\mu,\Sigma)d\theta\,.$
(4.26)

The quality of the Laplace approximation is determined by the overlap of the
typical set of the target distribution with that of the Gaussian
approximation. The greater the overlap, the more accurate the approximation
will be.

#### 21 Monte Carlo estimators

Another type of estimators is stochastic estimators, and in particular Monte
Carlo estimators [21] that rely on a collection of independent points or
samples,

$\\{\theta_{1},\dots,\theta_{n}\\}\in\Theta\,,$ (4.27)

from the distribution $p(\theta)$, such that the ensemble average of a
function $f(\theta)$,

$\hat{f}_{n}^{MC}=\frac{1}{n}\sum_{i=1}^{n}f(\theta_{n})\,,$ (4.28)

asymptotically converges to the corresponding expectation value

$\lim_{n\to\infty}\hat{f}_{n}^{MC}=\mathbb{E}_{p}[f(\theta)]\,.$ (4.29)

The asymptotic result of equation 4.29 is not particularly useful as a
computational algorithm will never be able to produce infinite samples.
Fortunately, the behaviour of Monte Carlo estimators can be quantified even
for finite samples.

For any square–integrable (i.e. both $\mathbb{E}_{p}[f]$ and
$\mathbb{E}_{p}[f^{2}]$ exist and are finite) real–valued function
$f(\theta)$, the Monte Carlo estimator satisfies the central limit theorem,

$\frac{\hat{f}_{n}^{MC}-\mathbb{E}_{p}[f(\theta)]}{\textrm{MC--
SE}_{n}[f]}\sim\mathcal{N}(0,1)\,,$ (4.30)

where $\textrm{MC--SE}_{n}[f]$ is the Monte Carlo Standard Error defined as,

$\textrm{MC--SE}_{n}[f]=\sqrt{\frac{\textrm{Var}_{p}[f]}{n}}\,.$ (4.31)

This means that we can estimate the expected number of samples that is
required to reach a certain level of precision for our estimates.

Figure 4.11: Illustration of the typical set including samples generated using
exact Monte Carlo sampling.

Another interesting property of Monte Carlo estimators is that their
precision, as quantified by the Monte Carlo Standard Error of equation 4.31,
does not depend on the dimensionality of the problem but relies only on the
number $n$ of samples instead. This means that Monte Carlo estimators can be
applied even in high–dimensional problems. This insensitivity to the curse of
dimensionality is directly related to the fact that the Monte Carlo samples
are already distributed in the typical set as shown in Figure 4.11. As we will
discover shortly, once we discuss more advanced methods, the difficult part is
to get the samples to the typical set in the first place.

We can also think of the Monte Carlo samples as a stochastic grid where the
computation is mostly focused in the regions of parameter space that
contribute to the computation of the expectation value. Starting from the
Monte Carlo estimator,

$\hat{f}_{n}^{MC}=\frac{1}{n}\sum_{i=1}^{n}f(\theta_{n})\,,$ (4.32)

and manipulate it into the quadrature form

$\hat{f}_{n}^{MC}=\sum_{i=1}^{n}f(\theta_{i})p(\theta_{i})\frac{1}{np(\theta_{i})}\,,$
(4.33)

where $\Delta\theta_{i}=1/np(\theta_{i})$ is the effective volume of each
sample.

Monte Carlo estimators are very powerful methods assuming that one can
generate independent samples from the target distribution. However, in most
interesting and realistic cases, this is not feasible. In that case, one has
to rely to alternative methods.

#### 22 Importance sampling

One alternative method to exact Monte Carlo sampling, that does not rely on
exact samples from the target distribution but instead requires an auxiliary
distribution is importance sampling. Importance sampling estimators use
samples from the auxiliary distribution and correct for any deviation from the
typical set of the target distribution using weighting factors. Although [105]
is typically credited with introducing importance sampling to statistics,
there are references to it in statistical physics as early as 1949 [71, 94].

In order to derive the importance weights necessary for the computation of the
expectation values we start with the definition of the expectation value and
do some re–arrangements,

$\begin{split}\mathbb{E}_{p}[f(\theta)]&=\int_{\Theta}f(\theta)p(\theta)d\theta\\\
&=\int_{\Theta}f(\theta)\frac{p(\theta)}{q(\theta)}q(\theta)d\theta\\\
&=\mathbb{E}_{q}\left[f(\theta)\times\frac{p(\theta)}{q(\theta)}\right]\,.\end{split}$
(4.34)

We can now estimate the expectation value,

$\hat{f}_{n}^{IS}=\frac{1}{n}\sum_{i=1}^{n}w(\tilde{\theta}_{i})f(\tilde{\theta}_{i})\,,$
(4.35)

using samples from the auxiliary distribution,

$\\{\tilde{\theta}_{1},\dots,\tilde{\theta}_{n}\\}\sim q(\tilde{\theta})\,,$
(4.36)

and importance weights given by,

$w(\tilde{\theta})=\frac{p(\tilde{\theta})}{q(\tilde{\theta})}\,.$ (4.37)

For any square–integrable real–valued function $f(\theta)$ , the importance
sampling estimator satisfies the central limit theorem,

$\frac{\hat{f}_{n}^{IS}-\mathbb{E}_{p}[f(\theta)]}{\textrm{IS--
SE}_{n}[f]}\sim\mathcal{N}(0,1)\,,$ (4.38)

where $\textrm{MC--SE}_{n}[f]$ is the Importance Sampling Standard Error
defined as,

$\textrm{IS--SE}_{n}[f]=\sqrt{\frac{\textrm{Var}_{q}[wf]}{n}}\,.$ (4.39)

By comparing the expression 4.39 for IS–SE to the respective expression 4.31
for the Monte Carlo Standard Error we can define the Effective Sample Size
(ESS),

$\textrm{ESS}_{n}[f]=\frac{\textrm{Var}_{q}[wf]}{\textrm{Var}_{p}[f]}n\,,$
(4.40)

as the effective number of exact samples that contain the same amount of
information as the $n$ samples and their importance weights.

It is important to mention here that if the target or auxiliary distribution
is known only up to a normalisation factor, for instance if the computation of
the normalisation constant is very costly, then the importance weights have to
be normalised such that,

$\sum_{i=1}^{n}w(\tilde{\theta}_{i})=1\,,$ (4.41)

for the aforementioned estimators to be valid.

Figure 4.12: Illustration of the typical set including samples generated using
an unsuitable importance density. In high dimensions the typical set
corresponds to a very thin shell and it is difficult to achieve sufficient
overlap between the typical set of the auxiliary and target distribution.
Here, this is depicted by samples that do not reside in the typical set of the
target and will thus have low importance weights.

The quality of the importance sampling estimator is determined by the amount
of overlap between the auxiliary and target distribution. Samples from the
auxiliary distribution residing in regions of high overlap will receive large
importance weights and those residing in regions of little or no overlap will
receive small importance weights.

As we showed in previous sections, in low dimensions the typical set is broad
so we should expect that the construction of importance sampling estimators
for low dimensional cases to be a feasible procedure. In higher dimensions
however, the typical shell is very thin thus complicating the choice of
auxiliary distributions with significant overlap with the target distribution.

It is common in practice to assume that an auxiliary distribution with broader
density tails than the target distribution would be sufficient to construct an
importance sampling estimator. Although the reasoning of this idea is
appealing, it can however be misleading as it does not extend to higher
dimensions in which the typical set becomes the central object of interest and
not the probability density.

It is useful to define a measure of the quality of an importance sampling
estimator. A straightforward choice would be to define the importance sampling
effective sample size,

$N_{\textrm{eff}}=\frac{\left(\sum_{i=1}^{n}w(\tilde{\theta}_{i})\right)^{2}}{\sum_{i=1}^{n}w(\tilde{\theta}_{i})^{2}}\,.$
(4.42)

Equation 4.42 is just a heuristic diagnostic and it should not be confused
with equation 4.40.

#### 23 Markov chain Monte Carlo

Importance sampling estimators trade the ability to produce exact samples from
the target distribution with weighted samples from an auxiliary distribution.
On the other hand, Markov chain Monte Carlo (MCMC) estimators replace the
exact samples with correlated samples generated by a Markov chain [68, 21].

Therefore, the key idea in MCMC is to explore the typical set using a sequence
of local steps. Starting a point $\theta_{1}$ in parameter space $\Theta$, the
next point $\theta_{2}$ is chosen stochastically in the neighbourhood of
$\theta_{1}$. Then the process is repeated for the next point $\theta_{3}$ in
the neighbourhood of $\theta_{2}$ and so on. At the end, we have generated a
chain of $n$ samples that is Markov, meaning that each sample conditionally
depends only on the previous one,

$P(\theta_{n}|\theta_{1},\dots,\theta_{n-1})=P(\theta_{n}|\theta_{n-1})\,.$
(4.43)

More formally, the Markov chain can be generated by repeatedly sampling from a
conditional probability distribution on the product space
$\Theta\times\Theta$, known as Markov transition probability
$T(\theta^{\prime}|\theta)$. Given an initial point $\theta_{1}$, sampling
from the Markov transition probability $T(\theta^{\prime}|\theta_{1})$ returns
sample $\theta_{1}$. We can thus construct a sequence of transitions,

$\begin{split}\theta_{2}&\sim T(\theta_{2}|\theta_{1})\\\ \theta_{3}&\sim
T(\theta_{3}|\theta_{2})\\\ &\dots\\\ \theta_{n}&\sim
T(\theta_{n}|\theta_{n-1})\,,\end{split}$ (4.44)

that constitute the Markov chain
$\\{\theta_{1},\theta_{2},\dots,\theta_{n}\\}$. The samples of the Markov
chain are not independent, but they are correlated. The reason for this is
their sequential origin i.e. $\theta_{n}$ depends on $\theta_{n-1}$ which
depends on $\theta_{n-2}$ so even samples that are not right next to each
other in the Markov chain can be correlated.

Figure 4.13: Example of trace plot of a Markov chain for parameter $\theta$.
The chain reaches the stationary state after about $50$ iterations.

##### 23.1 Requirements of MCMC

In general, the possible values $\theta$ of the Markov chain are called the
states of the Markov chain [172, 117]. For a Markov chain Monte Carlo
estimator to generate samples from the target distribution, the Markov chain
must satisfy a couple of conditions:

  1. 1.

First of all, the Markov chain has to leave the target distribution $p$
invariant or stationary,

$p(\theta^{\prime})=\int T(\theta^{\prime}|\theta)p(\theta)d\theta\,.$ (4.45)

This means that if we start from a state $\theta$ of $p$, the next state
$\theta^{\prime}$ is also a state of $p$. In practice, a sufficient but not
necessary condition is detailed balance, which requires that each transition
$\theta\rightarrow\theta^{\prime}$ is reversible. More formally, for any pair
of states $\theta$ and $\theta^{\prime}$ the following relation must hold,

$T(\theta^{\prime}|\theta)p(\theta)=T(\theta|\theta^{\prime})p(\theta^{\prime})\,,$
(4.46)

meaning that the probability of being at state $\theta$ and transitioning to
state $\theta^{\prime}$ is equal to the probability of being at state
$\theta^{\prime}$ and transitioning to state $\theta$.

  2. 2.

Furthermore, we need to make sure that the stationary distribution is unique
and that the distribution of states is able to converge to it regardless the
starting point $\theta_{1}$. In other words, we need to make sure that the
stationary distribution is also the limiting distribution. This requires two
properties, irreducibility, that is the ability to visit any state $\theta$
for which $p(\theta)>0$ in a finite number of steps, and aperiodicity, meaning
that no states are only accessible at certain regularly spaced times. These
two properties combined, when met, render the Markov chain ergodic.

##### 23.2 Expected behaviour

When all of the aforementioned conditions are obeyed, the Markov chain samples
from the target distribution. The behaviour of the Markov chain, in terms of
the computed expectation values, passes through four stages that characterise
its normal behaviour [15].

Figure 4.14: Initial stage of exploration – no exploration has taken place.

Figure 4.15: Burn–in stage of exploration – the absolute difference between
the estimate of $f$ and its expectation value slowly decreases as the chain
approaches the typical set.

The first stage, shown in Figure 4.14 consists of the initialisation of the
Markov chain. Often we do not know where the typical set resides and thus we
set the first state of the Markov chain to some arbitrary point in parameter
space. During the second stage, the Markov chain moves towards the typical set
as shown in Figure 4.15. At the same time the absolute difference of the
estimated value $\hat{f}$ from the expectation value $\mathbb{E}[f]$ slowly
decreases. In the third state shown in Figure 4.16, the Markov chain starts to
explore the typical set. The absolute difference between the estimate of $f$
and its expectation value decreases very rapidly. Finally, in the fourth stage
shown in Figure 4.17 the Markov chain wanders inside the typical set and the
standard error of the estimate asymptotically decreases as prescribed by the
central limit theorem.

Figure 4.16: Initial convergence phase – the absolute difference between the
estimate of $f$ and the expectation value decreases rapidly as the chain
approaches explores the typical set for the first time.

Figure 4.17: Stationary or equilibrium phase – the absolute difference between
the estimate of $f$ and the expectation value has reached its minimum value as
the chain samples fully populate the typical set.

##### 23.3 Central limit theorem of MCMC

Markov chain Monte Carlo estimators are particularly useful for many analyses
since they obey a central limit theorem that allows us to quantify their
precision [103, 67, 172]. In particular, given a square–integrable real–valued
function $f(\theta)$ and a long enough Markov chain, the following is true,

$\frac{\hat{f}_{n}^{\textrm{MCMC}}-\mathbb{E}_{p}[f(\theta)]}{\textrm{MCMC--
SE}_{n}[f]}\sim\mathcal{N}(0,1)\,,$ (4.47)

where $\textrm{MCMC--SE}_{n}[f]$ is the Markov chain Monte Carlo Standard
Error given by,

$\textrm{MCMC--
SE}_{n}[f]=\sqrt{\frac{\textrm{Var}_{p}[f]}{\textrm{ESS}_{n}[f]}}\,.$ (4.48)

Comparing $\textrm{MCMC--SE}_{n}[f]$ with the standard error $\textrm{MC--
SE}_{n}[f]$ of the exact Monte Carlo estimator given by equation 4.31, one
immediately notices that the number of samples $n$ has been replaced by the
term $\textrm{ESS}_{n}[f]$. This term, called the Effective Sample Size
accounts for the loss of information due to the correlation between samples
due to the Markov property of the chain. The effective sample size is given
by,

$\textrm{ESS}_{n}[f]=\frac{n}{\tau[f]}\,,$ (4.49)

where $\tau[f]$ is the relaxation or autocorrelation time of the Markov chain.
Less formally, $\tau[f]$ describes the number of steps required for the Markov
chain to “forget” where it started, meaning that only one out of $\tau[f]$ is
actually independent. A method for estimating the autocorrelation time of a
Markov chain will be discussed in the next subsection.

##### 23.4 Autocorrelation

The autocorrelation of the Markov chains is a necessary evil of MCMC methods
and must be properly understood before making any inference [67]. Figure 4.18
shows two Markov chains with different levels of autocorrelation. In the
weakly correlated chain, large jumps take place from one iteration to the
next. On the other hand, the strongly correlated chain is characterised by
very short jumps and more rigid trajectories.

Figure 4.18: Markov chains with different degrees of autocorrelation.

In order to quantify and measure the degree of autocorrelation of a Markov
chain we need to compare the states of the chain after fixed number of
iterations called lags. Given an arbitrary function $f(\theta)$ of the states,

$\mu_{f}=\mathbb{E}_{p}[f]\,,$ (4.50)

is the mean value of the function expressed as the expectation value over the
stationary distribution $p(\theta)$. The value $f(\theta_{i})-\mu_{f}$ then
quantifies the deviation of the $i$–th state of the chain from the mean value
$\mu_{f}$. The expectation value of the product of two such deviations defines
the autocovariance of the chain,

$c_{ij}=\mathbb{E}_{p}\left[(f_{i}-\mu_{f})(f_{j}-\mu_{f})\right]\,,$ (4.51)

where $f_{i}=f(\theta_{i})$ and $f_{j}=f(\theta_{j})$. Once the Markov chain
has reached the stationary phase, the autocovariance will no longer depend on
the particular states, $\theta_{i}$ and $\theta_{j}$, that we are comparing
but on the number of iterations, called lag $\ell=j-i$, between them,

$c_{ij}=c_{i,i+\ell}=c_{\ell}\,.$ (4.52)

Finally, if we normalise the autocovariance by the variance,

$\mathrm{Var}_{p}[f]=\mathbb{E}_{p}\left[(f-\mu_{f})^{2}\right]\,,$ (4.53)

we get the lag–$\ell$ autocorrelation function,

$\rho_{\ell}[f]=\frac{\mathbb{E}_{p}\left[(f_{i+\ell}-\mu_{f})(f_{i}-\mu_{f})\right]}{\mathrm{Var}_{p}[f]}\,.$
(4.54)

The normalisation ensures that the maximum possible value of $\rho_{\ell}[f]$
is $+1$ for fully correlated states and $-1$ for the completely
anti–correlated states. The value of $0$ corresponds to uncorrelated samples.
The lag–$\ell$ autocorrelation is always unity as any state is perfectly
correlated with itself. Furthermore, the autocorrelation function depends only
on the absolute lag and is invariant under changes of the sign, that is,
$\rho_{\ell}[f]=\rho_{-\ell}[f]$. For this reason, only the non–negative part
of the autocorrelation function is often plotted.

Figure 4.19: Autocorrelation as a function of lag $\ell$ for a weakly and a
strongly correlated chain.

Figure 4.19 shows the autocorrelation function for the weakly and strongly
correlated chains of Figure 4.18. We notice that although both functions begin
at the value of $1$ for lag $\ell=0$, they approach the value of $0$ at
different rates. In particular, the autocorrelation function of the weakly
correlated chain goes to $0$ after only a few lags, whereas the one
corresponding to the strongly correlated chain takes much longer.

The asymptotic variance of an infinitely long chain is defined as,

$\lim_{n\to\infty}n\times\mathrm{Var}_{p}\left[\hat{f}_{n}^{\mathrm{MCMC}}\right]=\mathrm{Var}_{p}[f]\times\tau_{I}[f]\,,$
(4.55)

where,

$\tau_{I}[f]=\sum_{\ell=-\infty}^{\ell=+\infty}\rho_{\ell}[f]=1+2\times\sum_{\ell=1}^{\ell=+\infty}\rho_{\ell}[f]\,,$
(4.56)

is the integrated autocorrelation time and the last equality hold due to the
lag–sign invariance of the autocorrelation function. Equation 4.55 implies
that the standard error of MCMC is,

$\textrm{MCMC--
SE}_{n}[f]=\sqrt{\frac{\mathrm{Var}_{p}[f]}{\mathrm{ESS}_{n}[f]}}\,,$ (4.57)

where we have defined the effective sample size as,

$\mathrm{ESS}_{n}[f]=\frac{n}{1+2\times\sum_{\ell=1}^{\ell=+\infty}\rho_{\ell}[f]}\,.$
(4.58)

Estimating the integrated autocorrelation time, and thus the effective sample
size, in not trivial in practice. The autocorrelation function can be very
noisy in large lags, as shown in Figure 4.19. This means that the sum in
equation 4.56 needs to be truncated in practice in order avoid adding noise.

### 5 Simple MCMC methods

Not all those who wander are lost.
— J.R.R. Tolkien

During the first half of the twentieth century, research efforts were focused
on the task of understanding the equilibrium behaviour of thermodynamic
systems. Furthermore, it was well understood that this behaviour was described
by specific probability distributions (e.g. canonical distribution for a
system in constant temperature). The physicists of that time showed great
interest in methods that produced exact samples from such probability
distributions. Enrico Fermi, for instance, would exploit such methods to make
amazingly–quick predictions of experimental outcomes as early as 1930s [121].
During the next two decades, Stan Ulam and John von Neumann developed various
such algorithms which, collectively, were anointed with the name “Monte Carlo”
after the infamous casino.

After the war, Nicholas Metropolis lead the group in Los Alamos in applying
Monte Carlo methods to increasingly complex thermodynamic systems. As exact
sampling was possible only for a limited number of distributions, Metropolis,
along with Arianna Rosenbluth, Marshall Rosenbluth, Edward Teller, and Augusta
Teller introduced the so–called Metropolis algorithm that produced correlated
samples from a wider range of probability distributions [122]. Arianna
implemented the algorithm on the MANIAC computer [75] and thus she is
considered the first person in history to implement a MCMC method.

After decades of empirical success in physics and chemistry, the statistician
Hastings [80] generalised the method by realising that by introducing a small
modification he could allow for any proposal distribution, not just a limited
family of symmetric ones. The method is known today as Metropolis–Hastings.
Despite Hasting’s seminal contribution, it was not until 1984 that [63]
introduced the Gibbs sampler for the task of image reconstruction, for the
broader statistical community to realise the potential of MCMC methods for
parameter inference [57].

#### 24 Metropolis–Hastings

The key idea of the Metropolis–Hastings algorithm is to separate the Markov
transition probability into two steps, a proposal and an acceptance step.
During the proposal step, a new state $\theta^{\prime}$ is generated
conditional on the current state $\theta$,

$\theta^{\prime}\sim q(\theta^{\prime}|\theta)\,,$ (5.1)

by sampling from a proposal distribution $q(\theta^{\prime}|\theta)$. The aim
of this step is to produce a new state that is likely, but not necessary, to
reside in the typical set of the target distribution. The form of the
conditional proposal distribution can be chosen based on the particular target
distribution. As we will see shortly, many of the developments in the field of
MCMC focus explicitly on the choice of the proposal distribution.

Once the new state $\theta^{\prime}$ is generated, its validity (i.e. whether
or not it belongs to the typical set) is assessed in the acceptance step. In
particular, the new state $\theta^{\prime}$ is accepted with probability,

$\alpha(\theta^{\prime},\theta)=\min\left(1,\frac{p(\theta^{\prime})q(\theta|\theta^{\prime})}{p(\theta)q(\theta^{\prime}|\theta)}\right)\,.$
(5.2)

Equation 5.2 is often called the Metropolis acceptance probability. If
accepted, the new state $\theta^{\prime}$ is added to the Markov chain and the
process is repeated with that as the current state (i.e.
$\theta\leftarrow\theta^{\prime})$. On the other hand, if the state
$\theta^{\prime}$ is rejected, the current state $\theta$ is added (i.e.
repeated) on the chain. This acceptance/rejection procedure based on equation
5.2 is often referred to as the Metropolis acceptance criterion.

It is important to mention here that the Metropolis acceptance criterion can
be evaluated even if we are only able to compute $p(\theta)$ up to a
normalisation constant, as any such factor would cancel out in the ratio
$p(\theta^{\prime})/p(\theta)$ that appears in equation 5.2. This is a very
important feature of the algorithm and one of the reasons for its widespread
success. In practice, it is very difficult to know the exact value for the
model evidence $\mathcal{Z}=p(d)$ that acts as the normalisation factor for
the posterior distribution $p(\theta|d)$ that might be the target
distribution.

It is straightforward to show that the Metropolis–Hastings algorithm leaves
the target distribution $p(\theta)$ stationary by first proving that it
preserves detailed balance. The Markov transition probability is simply,

$T(\theta^{\prime}|\theta)=q(\theta^{\prime}|\theta)\alpha(\theta^{\prime},\theta)\,,$
(5.3)

that is, the probability of proposing the new state $\theta^{\prime}$ given
the old state, times the probability of accepting it. Therefore, the Markov
transition probability for the Metropolis–Hastings algorithm preserves
detailed balance,

$\begin{split}T(\theta^{\prime}|\theta)p(\theta)&=q(\theta^{\prime}|\theta)\alpha(\theta^{\prime},\theta)p(\theta)\\\
&=q(\theta^{\prime}|\theta)\min\left(1,\frac{p(\theta^{\prime})q(\theta|\theta^{\prime})}{p(\theta)q(\theta^{\prime}|\theta)}\right)p(\theta)\\\
&=\min\left(p(\theta)q(\theta^{\prime}|\theta),p(\theta^{\prime})q(\theta|\theta^{\prime})\right)\\\
&=q(\theta|\theta^{\prime})\alpha(\theta,\theta^{\prime})p(\theta^{\prime})\\\
&=T(\theta|\theta^{\prime})p(\theta^{\prime})\,,\end{split}$ (5.4)

Algorithm 2 Metropolis–Hastings

0: initial state $\theta_{1}$, (unnormalised) target density $f(\theta)\propto
p(\theta)$, proposal density $q(\theta^{\prime}|\theta)$, and number of
iterations $N$

0: Markov chain $\theta_{1},\theta_{2},\dots,\theta_{N}$ that has $p(\theta)$
as its equilibrium distribution

1: for $t=1$ to $N$ do

2: Draw new state from proposal distribution $\theta^{\prime}\sim
q(\theta^{\prime}|\theta_{t})$

3: Compute acceptance probability
$\alpha=\min\left(1,\frac{f(\theta^{\prime})q(\theta_{t}|\theta^{\prime})}{f(\theta_{t})q(\theta^{\prime}|\theta_{t})}\right)$

4: Draw uniform random number $u\sim\mathcal{U}(0,1)$

5: if $u<\alpha$ then

6: Accept proposed state and set $\theta_{t+1}\leftarrow\theta^{\prime}$

7: else

8: Reject proposed state and set current state as the next
$\theta_{t+1}\leftarrow\theta_{t}$

9: end if

10: end for

##### 24.1 Random–walk Metropolis

A very common, and simplifying in practice, choice for the proposal
distribution is the conditional normal distribution
$q(\theta^{\prime}|\theta)=\mathcal{N}(\theta^{\prime}|\theta,\Sigma)$ centred
around the current state $\theta$ with covariance matrix $\Sigma$ [122, 172].
The probability density has the usual Gaussian functional form,

$q(\theta^{\prime}|\theta)=\det(2\pi\Sigma)^{-\frac{1}{2}}\exp{\left[-\frac{1}{2}(\theta^{\prime}-\theta)^{T}\Sigma^{-1}(\theta^{\prime}-\theta)\right]}\,.$
(5.5)

The symmetry of this proposal distribution,

$q(\theta^{\prime}|\theta)=q(\theta|\theta^{\prime})\,,$ (5.6)

means that the Metropolis acceptance probability of equation 5.2 is simplified
and the $q$ terms drop out,

$\alpha(\theta^{\prime},\theta)=\min\left(1,\frac{p(\theta^{\prime})}{p(\theta)}\right)\,.$
(5.7)

##### 24.2 Independence Metropolis

Another simple choice of proposal distribution is to make it independent of
the current state $\theta$. For instance, one can choose a normal distribution
$q(\theta)=\mathcal{N}(\theta|\mu,\Sigma)$ with mean $\mu$ and covariance
matrix $\Sigma$, both of which must be known a priori and can not depend on
the current state. In this case, the Metropolis acceptance probability reduces
to,

$\alpha(\theta^{\prime},\theta)=\min\left(1,\frac{p(\theta^{\prime})q(\theta)}{p(\theta)q(\theta^{\prime})}\right)\,.$
(5.8)

One of the benefits of Independence Metropolis [80, 172], as this approach is
called, is that any states produced are independent samples from the target
distribution $p(\theta)$. However, it suffers from similar problems to
Importance sampling. Instead of vanishingly small importance weights, in
Independence Metropolis we might experience vanishingly small acceptance
probability when the overlap of the typical set of the proposal distribution
$q$ with the target $p$ is small. For this reason, the use of Independence
Metropolis is wise only when we have good reasons to believe that the proposal
distribution is sufficiently close to the target distribution or the
dimensionality is low.

##### 24.3 Metropolis–adjusted Langevin algorithm

As we have seen, the normal proposal distribution of Random–walk Metropolis
can utilise only global information about the target distribution (i.e. the
covariance matrix) in order to achieve efficient sampling. Although sufficient
in low to moderate dimensional problems, this strategy can become inefficient
as the number of parameters of the target distribution increases. In practice,
Random–walk Metropolis proposes new states indiscriminately along directions
of great covariance without taking into account the local structure of the
typical set. This results in either low acceptance probabilities or small
proposed steps being accepted as the typical set become thinner in high
dimensions.

One way to circumvent this effect and achieve better sampling performance is
to capitalise on the knowledge of the gradient of the target distribution in
order to “bias” the proposed states towards directions that are more likely to
lead to higher acceptance probability. Metropolis–adjusted Langevin algorithm
(MALA) [153] achieves this by using a conditional normal distribution,

$q(\theta^{\prime}|\theta)=\mathcal{N}\left(\theta+\tau\Sigma\nabla\log
p(\theta),2\tau\Sigma\right)\,,$ (5.9)

where its mean $\theta+\tau\Sigma\nabla\log p(\theta)$ is shifted from the
current state $\theta$ along the direction of the gradient of the logarithm of
the target distribution $\nabla\log p(\theta)$. If known, $\Sigma$ can be an
approximate covariance matrix that characterises the target distribution,
otherwise, a unit–diagonal matrix can be used. $\tau$ is the step size of the
method and determines the amount of shift of the proposal distribution. In the
limit that $\tau\rightarrow 0$, MALA reduces to RWM. The step size $\tau$ can
be modified in order to achieve the theoretically optimal acceptance
probability of $0.574$. Despite the fact that the aforementioned acceptance
rate has only been proven to be optimal for certain types of target
distributions [151], we expect that values in the range between $0.4$ and
$0.8$ would result in a high performance for most applications.

In terms of the typical set, we can think of the gradient of the log
probability as a guide that allows for better–informed proposals that are more
likely to belong to the typical shell.

##### 24.4 Adaptive Metropolis

Hyperparameters of the proposal distribution, such as the covariance matrix
$\Sigma$ of RWM or the step size $\tau$ of MALA, do not have to be chosen a
priori or based on preliminary MCMC runs but they can instead be adaptively
tuned during the run. [77] presented a prototype adaptive version of RWM in
which the proposal distribution is continuously adapted during the run using
all of the collected samples in order to estimate its covariance matrix. The
estimation of the covariance matrix is efficient as only incremental updates
are required using simple recursive formulas.

In order to achieve this kind of proposal adaptation in practice we need to
abandon the Markov property of the chain. In general, this is not a problem as
there is nothing special about the Markov property apart from its simplicity.
However, continuous tuning of the proposal distribution during the run
requires the adaptation to be diminishing, with very specific characteristics,
in order to preserve the ergodicity of the method [21].

One of the most commonly used algorithms for diminishing adaptation is the
stochastic approximation algorithm of [148]. Suppose that we have a function
$f(\lambda_{i})=f_{i}$ , which encodes some aspect of the behaviour of the
$i$–th state of chain (e.g. the acceptance probability) as a function of some
tunable property $\lambda$ (e.g. the proposal scale), that has expectation
value,

$\mathbb{E}[f(\lambda)]=\frac{1}{n}\sum_{i=1}^{n}f(\lambda_{i})\,.$ (5.10)

The solution to the equation $\mathbb{E}[f(\lambda)]=f^{*}$ can be found
iteratively, using the recursive formula,

$\lambda_{i+1}=\lambda_{i}-\gamma_{i}\left(f_{i}-f^{*}\right)\,,$ (5.11)

assuming that $f$ is a non–decreasing function of $\lambda$ that is uniformly
bounded [2]. The parameter $\gamma_{i}$ determines the learning rate or the
rate of convergence of the approximation and has to obey two conditions,

$\sum_{i=1}^{n}\gamma_{i}=\infty\,,\quad\sum_{i=1}^{n}\gamma_{i}^{2}<\infty\,.$
(5.12)

The former condition ensures that any point $\theta$ can eventually be
reached, and the latter condition ensures that the fluctuations introduced by
new iterations is contained and does not prevent convergence. A commonly used
schedule for the learning rate that satisfies the above conditions has the
form $\gamma_{i}=i^{-\kappa}$ for $\kappa\in(0.5,1]$.

Algorithm 3 Adaptive Metropolis

0: initial state $\theta_{1}$, (unnormalised) target density $f(\theta)\propto
p(\theta)$, the target acceptance rate $\alpha^{*}$, learning rate schedule
(e.g. $g_{t}=1/t$), and number of iterations $N$

0: Markov chain $\theta_{1},\theta_{2},\dots,\theta_{N}$ that has $p(\theta)$
as its equilibrium distribution

1: Initialise $\mu_{1}=0$, $\Sigma_{1}=1$, $\log\lambda_{1}=0$

2: for $t=1$ to $N$ do

3: Draw new state from proposal distribution
$\theta^{\prime}\sim\mathcal{N}(\theta^{\prime}|\theta_{t},\lambda_{t}\Sigma_{t})$

4: Compute acceptance probability
$\alpha_{t}=\min\left(1,f(\theta^{\prime})/f(\theta_{t})\right)$

5: Draw uniform random number $u\sim\mathcal{U}(0,1)$

6: if $u<\alpha_{t}$ then

7: Accept proposed state and set $\theta_{t+1}\leftarrow\theta^{\prime}$

8: else

9: Reject proposed state and set current state as the next
$\theta_{t+1}\leftarrow\theta_{t}$

10: end if

11: Update mean estimate
$\mu_{t+1}\leftarrow\mu_{t}-\gamma_{t}(\mu_{t}-\theta_{t+1})$

12: Update covariance estimate
$\Sigma_{t+1}\leftarrow\Sigma_{t}-\gamma_{t}\left[\Sigma_{t}-(\mu_{t}-\theta_{t+1})(\mu_{t}-\theta_{t+1})^{T}\right]$

13: Update proposal scale estimate
$\log\lambda_{t+1}\leftarrow\log\lambda_{t}-\gamma_{t}(\alpha_{t}-\alpha^{*})$

14: end for

Let us now go through an example of developing an adaptive version of the
commonly used RWM, in which we tune the covariance matrix $\Sigma$ of the
Gaussian proposal distribution,
$\theta^{\prime}\sim\mathcal{N}(\theta^{\prime}|\theta,\lambda\Sigma)$, using
the following diminishing adaptation scheme,

$\begin{split}\mu_{i+1}&=\mu_{i}-\gamma_{i}\left(\mu_{i}-\theta_{i+1}\right)\,,\\\
\Sigma_{i+1}&=\Sigma_{i}-\gamma_{i}\left[\Sigma_{i}-(\mu_{i}-\theta_{i+1})(\mu_{i}-\theta_{i+1})^{T}\right]\,,\end{split}$
(5.13)

where the $\mu$ is the mean value used for the estimation of the covariance
$\Sigma$, and $\gamma_{i}=1/i$ is the learning rate. At the same time we can
also tune the magnitude of the proposal scale, $\lambda_{i}$, by attempting to
match the acceptance probability $\alpha$ to the theoretically optimal value
of $\alpha^{*}=0.234$,

$\log\lambda_{i+1}=\log\lambda_{i}-\gamma_{i}\left(\alpha_{i}-\alpha^{*}\right)\,.$
(5.14)

Understanding equation 5.14 is straightforward, if the observed acceptance
rate is greater than the target (i.e. $\alpha_{i}<\alpha^{*}$) then the
logarithm of the magnitude of the proposal scale $\log\lambda$ is reduced and
vice versa. The Adaptive Metropolis method presented in this paragraph
constitutes a generalisation of the method presented by [77]. With the
inclusion of the adaptation of the proposal scale using equation 5.14 the
algorithm resembles that of [2].

#### 25 Gibbs sampling

Another very popular Markov chain Monte Carlo method, to which we partly owe
the widespread use of Bayesian inference today, is Gibbs sampling. Initially
known as the heat bath algorithm in the statistical physics literature, the
Gibbs sampler enjoyed great success in the statistical community following the
seminal paper by [63] that demonstrated its benefits for analysing Gibbs
distributions on lattices in the context of image processing.

##### 25.1 Gibbs sampler

Gibbs sampler attempts to overcome the curse of dimensionality using
conditioning [26]. In particular, assuming that exact sampling from the
conditional distributions of the target distribution $p(\theta)$ is possible,
we can generate samples from the target distribution by sequentially sampling
from its full set of conditionals instead. Given an initial state
$\theta=(\theta_{1},\dots,\theta_{D})$, the next state in the Markov chain can
be generated as,

$\begin{split}\theta_{1}^{\prime}&\sim
p(\theta_{1}|\theta_{2},\dots,\theta_{D})\\\ &\vdots\\\
\theta_{k}^{\prime}&\sim
p(\theta_{k}^{\prime}|\theta_{1},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots\theta_{D})\\\
&\vdots\\\ \theta_{D}^{\prime}&\sim
p(\theta_{D}^{\prime}|\theta_{1},\dots,\dots\theta_{D-1}^{\prime})\,.\end{split}$
(5.15)

The current state $\theta$ is then replaced by the new state
$\theta^{\prime}=(\theta_{1}^{\prime},\dots,\theta_{D}^{\prime})$ and the
process is repeated until enough states are collected in the Markov chain. The
order of the state updates of equation 5.15 can be either fixed (with a
possible reversal after every iteration), as shown above, or randomised to
ensure detailed balance.

Algorithm 4 Gibbs sampler

0: initial state
$\theta^{(1)}=(\theta_{1}^{(0)},\theta_{2}^{(0)},\dots,\theta_{D}^{(0)})$, all
the conditional distributions
$p(\theta_{k}|\theta_{1},\dots,\theta_{k-1},\theta_{k+1},\dots,\theta_{D})$ of
target $p(\theta)$, and number of iterations $N$

0: Markov chain $\theta^{(1)},\theta^{(2)},\dots,\theta^{(N)}$ that has
$p(\theta)$ as its equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $D$ do

3: Draw from conditional $\theta_{k}^{\prime}\sim
p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})$

4: end for

5: Set
$\theta^{(t+1)}\leftarrow\theta^{\prime}=(\theta_{1}^{\prime},\dots,\theta_{D}^{\prime})$

6: end for

##### 25.2 Metropolis–within–Gibbs sampler

The Gibbs sampler relies on our ability to produce samples from each one of
the conditional distributions. This however is not always feasible as some of
the components of the full conditional set might not admit an exact sampling
solution. Instead of abandoning Gibbs sampler altogether, [124, 125] suggested
the use of a compromise between the Gibbs sampler and the Metropolis–Hastings
algorithm.

The key idea behind the Metropolis–within–Gibbs sampler is to use Gibbs
sampling for as many of the conditional distributions as possible in order to
produce exact samples, and rely on correlated samples generated using
Metropolis–Hastings for any conditional distributions that exact sampling is
not possible.

Suppose that we have a partial state
$(\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k},\dots,\theta_{D}$)
and we have difficulty generating exact samples from the conditional
distribution $\theta_{k}^{\prime}\sim
p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})$.
We can then treat
$p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})$
as the target distribution for a Metropolis–Hastings estimator as follows, in
order to proceed with the computation,

  1. 1.

First, we have to propose a new sample $\theta_{k}^{*}\sim
q(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})$
from an arbitrary proposal distribution,

  2. 2.

Then, compute the Metropolis acceptance probability

$\begin{split}\alpha(\theta_{k}^{*},\theta_{k})=\min\bigg{(}1,&\frac{p(\theta_{k}^{*}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})}{p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})}\\\
\times&\frac{q(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})}{q(\theta_{k}^{*}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1},\dots,\theta_{D})}\bigg{)}\,,\end{split}$
(5.16)

  3. 3.

Finally, accept the new state $\theta_{k}^{\prime}\leftarrow\theta_{k}^{*}$
with probability $\alpha(\theta_{k}^{*},\theta_{k})$, otherwise reject and
keep the previous state $\theta_{k}^{\prime}\leftarrow\theta_{k}$ .

Using the algorithm presented above we can replace exact sampling from
conditional distributions where it is not feasible with Metropolis–Hastings
estimates.

Algorithm 5 Metropolis–within–Gibbs sampler

0: initial state
$\theta^{(1)}=(\theta_{1}^{(0)},\theta_{2}^{(0)},\dots,\theta_{D}^{(0)})$,
proposal distributions
$q(\theta_{k}|\theta_{1},\dots,\theta_{k-1},\theta_{k+1},\dots,\theta_{D})$,
(unnormalised) target distribution $f(\theta)\propto p(\theta)$, the
conditional distributions
$p(\theta_{\ell}|\theta_{1},\dots,\theta_{\ell-1},\theta_{\ell+1},\dots,\theta_{D})$
of target $p(\theta)$ where $\ell\in L$ and $L$ the set of indices for which
exact sampling of the conditional is possible, and number of iterations $N$

0: Markov chain $\theta^{(1)},\theta^{(2)},\dots,\theta^{(N)}$ that has
$p(\theta)$ as its equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $D$ do

3: if $k\in L$ then

4: Draw from conditional $\theta_{k}^{\prime}\sim
p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})$

5: else

6: Draw from proposal $\theta_{k}^{*}\sim
q(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})$

7: Compute acceptance probability
$\begin{split}\alpha=\min\bigg{(}1,&\frac{p(\theta_{k}^{*}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})}{p(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})}\\\
\times&\frac{q(\theta_{k}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})}{q(\theta_{k}^{*}|\theta_{1}^{\prime},\dots,\theta_{k-1}^{\prime},\theta_{k+1}^{(t)},\dots,\theta_{D}^{(t)})}\bigg{)}\end{split}$

8: Draw uniform number $u\sim\mathcal{U}(0,1)$

9: if $u<\alpha$ then

10: accept new partial state and set
$\theta_{k}^{\prime}\leftarrow\theta_{k}^{*}$

11: else

12: reject new partial state and set
$\theta_{k}^{\prime}\leftarrow\theta_{k}^{(t)}$

13: end if

14: end if

15: end for

16: Set
$\theta^{(t+1)}\leftarrow\theta^{\prime}=(\theta_{1}^{\prime},\dots,\theta_{D}^{\prime})$

17: end for

### 6 Auxiliary variable MCMC methods

Natura non facit saltus.
— Gottfried Leibniz

Auxiliary variable MCMC methods rely on the introduction of one or more
additional variables in order to make sampling from the target distribution
more efficient.

#### 26 Simulated annealing

In metallurgy, annealing refers to the thermal process used to harden steel.
Initially, the metal is heated to a high temperature and then it is cooled
down slowly enough for the atoms to self–arrange in an ordered pattern that
corresponds to the minimum energy [24]. The slow rate of cooling ensures that
the energy of the system will reach its global minimum instead of getting
trapped in local minima.

Realising that the Metropolis–Hastings method can be used to simulate the
process of gradually cooling a solid towards a low–temperature equilibrium
state, [104] suggested that we should construct a sequence of Boltzmann
distributions,

$p_{i}(\theta)\propto e^{-\frac{H(\theta)}{T_{i}}}\,,$ (6.1)

for a series of temperatures $T_{1}>T_{2}>\dots>T_{m}$ and simulate from each
one in succession by performing a number of MCMC steps in each temperature
before moving on to the next. The result is a non–homogeneous Markov chain,
that is, a Markov chain with time–varying target density. Assuming that
$T_{1}$ is high enough and $T_{m}\approx 0$ we can find the global minimum of
the energy $E(\theta)$ by simulating the cooling process of a solid.

Figure 6.1: Illustration of the gradual annealing performed in the posterior
distribution. The prior distribution corresponds to $1/T\rightarrow 0$ and the
posterior is recovered as $1/T\rightarrow 1$.

Simulated annealing can be used for sampling too, not just for optimisation.
By stopping the cooling process earlier at $T_{m}=1$ we can sample from any
target distribution $p(\theta)$, not just Boltzmann distributions, simply
setting $H(\theta)=-\log p(\theta)$. Furthermore, for applications in which
the target distribution is the posterior distribution we can construct the
following sequence of densities,

$p_{i}(\theta)\propto p(\theta)e^{\frac{\log p(d|\theta)}{T_{i}}}\,,$ (6.2)

where $T_{1}>T_{2}>\dots>T_{m}=1$. In this case, for $T_{1}>>1$ we effectively
sample from the prior distribution,

$p_{i}(\theta)\propto p(\theta)\,,$ (6.3)

whereas in the limit that $T_{m}=1$ we acquire samples from the posterior,

$p_{i}(\theta)\propto p(\theta)p(d|\theta)\,.$ (6.4)

The number of MCMC steps to perform in each temperature, before moving on to
the next one, is arbitrary and different mixing criteria can be utilised (e.g.
Gelman–Rubin, autocorrelation thresholds, etc.). The benefit of using
simulated annealing for sampling is that by simulating multiple Markov chains,
possibly in parallel, through this sequence of densities, the risk of the
chains getting trapped in isolated modes of the posterior distribution is
minimised. This means that this approach can be used when the probability
distribution is strongly multimodal. Furthermore, if the number of temperature
levels is large enough and the spacing between them small enough, then the
Markov chain is approximately always in equilibrium, meaning that no, or
minor, burn–in is required to be discarded.

Algorithm 6 Simulated annealing

0: initial state $\theta_{1}^{(1)}$, temperature schedule
$T_{1}>T_{2}>\dots>T_{m}=1$, prior density $\pi(\theta)\equiv
p(\theta|\mathcal{M})$, likelihood function $\mathcal{L}(\theta)\equiv
p(d|\theta,\mathcal{M})$, and number of MCMC iterations $N$ per temperature

0: Multiple Markov chains $\theta^{(1)},\theta^{(2)},\dots,\theta^{(m)}$ with
the last one having $p(\theta)$ as its equilibrium distribution

1: for $i=1$ to $m$ do

2: Set annealed density
$p_{i}(\theta)\propto\pi(\theta)\mathcal{L}(\theta)^{1/T_{i}}$

3: Generate Markov chain $\theta_{1}^{(i)},\dots,\theta_{N}^{(i)}$ targeting
$p_{i}(\theta)$ (e.g. using Metropolis–Hastings)

4: Set last state as the first state for the next annealed density
$\theta_{1}^{(i+1)}\leftarrow\theta_{N}^{(i)}$

5: end for

#### 27 Slice sampling

Slice sampling is another MCMC method that relies on an auxiliary variable in
order to make sampling easier [13, 131, 133]. The method is based on the
realisation that sampling from the target distribution with density
$p(\theta)$ is equivalent to uniform sampling from the area or volume below
the curve or surface of $f(\theta)\propto p(\theta)$. This is equivalent to
the introduction of an auxiliary variable $\phi$, called height, such that the
joint distribution $p(\theta,\phi)$ is uniform over the region,

$U=\left\\{(\theta,\phi)\,:\,0<\phi<f(\theta)\right\\}\,.$ (6.5)

In other words, the joint distribution can be written as,

$p(\theta,\phi)=\begin{cases}1/\mathcal{Z}&\text{if }0<\phi<f(\theta)\,,\\\
0&\text{otherwise}\,,\end{cases}$ (6.6)

where

$\mathcal{Z}=\int f(\theta)d\theta\,.$ (6.7)

To sample from the target distribution $p(\theta)$ we first sample uniformly
from $p(\theta,\phi)$ and then marginalise over $\phi$ by dropping the
$\phi$–value of each sample and keeping the $\theta$–value. The proof that
this results in the marginal density for $\theta$ is straightforward,

$p(\theta)=\int
p(\theta,\phi)d\phi=\int_{0}^{f(\theta)}\frac{1}{\mathcal{Z}}d\phi=\frac{f(\theta)}{\mathcal{Z}}\,.$
(6.8)

Generating independent samples from the uniform joint density $p(\theta,\phi)$
is rarely possible in practice. Instead, one might prefer to construct a
Markov chain that leaves the distribution $p(\theta,\phi)$ invariant. One such
option is to use Gibbs sampling, that is, to sample alternately from the
conditional distribution $p(\phi|\theta)$, which is uniform over the interval
$(0,f(\theta))$, and then from the conditional distribution $p(\theta|\phi)$,
which is uniform over the region,

$S=\left\\{\theta\,:\,\phi<f(\theta)\right\\}\,,$ (6.9)

called the slice. Applying this procedure repeatedly will produce a Markov
chain that has the joint distribution $p(\theta,\phi)$ as its stationary
distribution.

Sampling uniformly from the aforementioned slice is not trivial either.
However, the fact that the conditional density $p(\theta|\phi)$ is uniform
allows us to construct procedures to sample from it which would otherwise
would not have worked. [133] proposed the following sequence of steps for
univariate probability distributions,

  1. 1.

Uniformly sample a real value $\phi$ in the interval $(0,f(\theta_{0}))$,
therefore defining the horizontal slice
$S=\left\\{\theta\,:\,\phi<f(\theta_{0})\right\\}$ that always includes
$\theta_{0}$,

  2. 2.

Find an interval $I=(L,R)$ around $\theta_{0}$ along the slice that contains
all, or much of, the slice,

  3. 3.

Sample a new value $\theta_{1}$ from the part of the slice within the
interval, that is, from $I\cap S$.

It is important to mention that as we often work with $g(\theta)=\log
f(\theta)$, to avoid numerical issues, one can use the variable
$\psi=\log(\phi)=g(\theta_{0})-e$, where $e$ is exponentially distributed with
mean one, to define the slice as $S=\\{\theta\,:\,\psi<g(\theta_{0})\\}$.

Figure 6.2: Illustration of the stepping–out and shrinking procedures used in
slice sampling. Given an initial state $\theta$ in the Markov chain, an
auxiliary variable $\phi$ is sampled corresponding to the height thus defining
the extended state $(\theta,\phi)$ shown here as a blue point. An interval of
a certain width is placed uniformly around the current point $(\theta,\phi)$
and expanded in steps of size equal to the initial width until both of its
ends, $L$ and $R$, are outside the graph. A new state, shown in red, is then
proposed uniformly along the interval $(L,R)$. Since the proposed state lies
above the graph of $f(\theta)$ (i.e. not in the slice shown as a continuous
line) it is rejected. A new state, shown in green, is the proposed uniformly
between the rejected state and $R$. Since the proposed state is below the
graph, and thus in the slice, it is accepted and added to the Markov chain.
The whole process is then repeated.

The first step in the above procedure is trivial, yet steps two and three
require more serious consideration. [133] suggested to use the so–called
stepping–out and shrinking procedures for those steps respectively.
Stepping–out works by uniformly positioning an interval of width $w$ around
$\theta_{0}$ such that it includes $\theta_{0}$, and then expanding the
interval in steps of size $w$ until both its ends $(L,R)$ are outside the
slice $S$. This effectively constructs the interval $I=(L,R)$. It is worth
noting that the algorithm is valid even if only a pre–specified number of
expansions take place and the interval ends up not covering the entirety of
the slice $S$. The shrinking procedure that follows functions by uniformly
sampling points in the interval $I=(L,R)$ until one of them lies in the slice
$S$. Every time a point is rejected, being outside of the slice, the interval
$I$ shrinks such that the rejected point now defines one of its two
boundaries, determined by whether the rejected point lies left or right of
$\theta_{0}$.

The fact that the three–step procedure presented so far describes a slice
sampling update from a univariate probability distribution $p(\theta)$ does
not prohibit its use in multivariate cases. In particular, there are many ways
the aforementioned recipe can be generalised and used in target distribution
with more than one parameter. Perhaps the simplest one is to apply this
univariate scheme along each coordinate axis in turn, updating one parameter
at a time. This corresponds to a Metropolis–within–Gibbs scheme. Another
option is to apply 1–D updates in random directions. This is more general than
the previous one, and there is freedom to choose the distribution of the
random directions. The directions can be drawn from a multivariate zero–mean
normal distribution with unit–diagonal covariance matrix or a more appropriate
non–diagonal covariance matrix that encodes some of the correlations of the
parameters of the target distribution. Such a covariance can be configured a
priori, estimated during a short preliminary run from samples from the target,
or adaptively tuned using an appropriate algorithm for diminishing adaptation.

One of the great benefits of slice sampling is the fact that it has a single
hyper–parameter, the initial width $w$ of the interval $I$. Furthermore, the
value of $w$ is adapted continuously by the stepping–out and shrinking
procedures. This sort of local adaptation is absent from many MCMC that assume
a global proposal scale. Another characteristic of slice sampling is the lack
of rejected samples in the Markov chain. Unlike methods that include a
Metropolis acceptance criterion, slice sampling always moves to a new state in
every iteration.

Algorithm 7 Slice sampling

0: initial state $\theta_{1}$, (unnormalised) target density $f(\theta)\propto
p(\theta)$, number of maximum expansions $m$, and number of iterations $N$

0: Markov chain $\theta_{1},\dots,\theta_{N}$ that has $p(\theta)$ as its
equilibrium distribution

1: Draw “height” auxiliary variable $\phi\sim\mathcal{U}(0,f(\theta_{1}))$

2: for $t=1$ to $N$ do

3: Draw left bound of the $I$ interval
$L\sim\mathcal{U}(\theta_{t}-w,\theta_{t})$

4: Set right bound $R\leftarrow L+w$

5: Draw uniform variable $u\sim\mathcal{U}(0,1)$

6: Set maximum number of interval expansions to the left
$J\leftarrow\text{Floor}(m\times u)$

7: Set maximum number of interval expansions to the right $K\leftarrow(m-1)-J$

8: while $J>0$ and $\phi<f(L)$ do

9: Expand left boundary $L\leftarrow L-w$ in steps of $w$

10: Reduce count $J\leftarrow J-1$ by one

11: end while

12: while $K>0$ and $\phi<f(R)$ do

13: Expand right boundary $L\leftarrow R+w$ in steps of $w$

14: Reduce count $K\leftarrow K-1$ by one

15: end while

16: repeat

17: Draw state within interval $\theta^{\prime}\sim\mathcal{U}(L,R)$

18: if $\theta^{\prime}<\theta_{t}$ then

19: Contract interval $L\leftarrow\theta^{\prime}$

20: else

21: Contract interval $R\leftarrow\theta^{\prime}$

22: end if

23: until $\phi<f(\theta^{\prime})$

24: Accept new state $\theta_{t+1}\leftarrow\theta^{\prime}$

25: end for

#### 28 Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) introduces a momentum auxiliary variable and
uses the gradient of the target probability density to efficiently explore the
typical set. HMC turns the problem of sampling from the target distribution
into the approximate simulation of Hamiltonian dynamics with a subsequent
Metropolis correction step [135]. In the statistical physics literature HMC
was suggested as a method of efficiently simulating states from a physical
system [39], which was then employed to statistical inference problems [128,
129, 130, 116].

Figure 6.3: Illustration of Hamiltonian trajectories in parameter space. The
black points correspond to the accepted states.

##### 28.1 Auxiliary momentum variable

HMC introduces an auxiliary variable $\rho$ and samples from the joint
probability density,

$p(\rho,\theta)=p(\rho|\theta)p(\theta)\,.$ (6.10)

In most applications of HMC, the momentum variable $\rho$ chosen to be
Gaussian–distributed,

$\rho\sim\mathcal{N}(0,M)\,,$ (6.11)

and its probability density function to be independent of the state variable
$\theta$ (i.e. $p(\rho|\theta)=p(\theta)$). $M$ is the symmetric, positive
definite mass matrix that has the role of the Euclidean metric, that is to
define the relative length scales between parameters. In practice, $M$ can be
chosen to be $\Sigma^{-1}$, meaning the inverse of the sample covariance
matrix that characterises the target distribution assuming that it is known or
easy to estimate.

##### 28.2 The Hamiltonian

HMC treats sampling from the joint distribution $p(\rho,\theta)$ as a problem
of solving the Hamiltonian dynamics given the Hamiltonian,

$\begin{split}\mathcal{H}(\rho,\theta)&=-\log p(\rho,\theta)\\\ &=-\log
p(\rho|\theta)-\log p(\theta)\\\ &=T(\rho|\theta)+V(\theta)\,,\end{split}$
(6.12)

where

$T(\rho|\theta)=-\log p(\rho|\theta)\,,$ (6.13)

is the kinetic energy, and,

$V(\theta)=-\log p(\theta)\,,$ (6.14)

is the potential energy.

##### 28.3 Hamilton’s equations

The dynamics of a system (i.e. its evolution in time) that is characterised by
the Hamiltonian of equation 6.12 are given by solving Hamilton’s equations,

$\begin{split}\frac{d\theta}{dt}&=\frac{\partial\mathcal{H}}{\partial\rho}=\frac{\partial
T}{\partial\rho}\,,\\\
\frac{d\rho}{dt}&=-\frac{\partial\mathcal{H}}{\partial\theta}=-\frac{\partial
T}{\partial\theta}-\frac{\partial V}{\partial\theta}\,,\end{split}$ (6.15)

or, in the case that the momentum variable $\rho$ is independent of the state
variable $\theta$, that is $p(\rho|\theta)=p(\rho)$,

$\begin{split}\frac{d\theta}{dt}&=\frac{\partial\mathcal{H}}{\partial\rho}=\frac{\partial
T}{\partial\rho}\,,\\\
\frac{d\rho}{dt}&=-\frac{\partial\mathcal{H}}{\partial\theta}=-\frac{\partial
V}{\partial\theta}\,.\end{split}$ (6.16)

Therefore, given an initial state $(\rho,\theta)$, the system’s evolution in
time is completely determined by equations 6.16.

##### 28.4 Leapfrog integration

Solving Hamilton’s equations analytically is only feasible for very simple
systems that correspond to simple target probability distributions. In
practice, however, we aim to solve equations 6.16 for systems of arbitrary
complexity. To this end, we turn to numerical methods for integrating this
system of differential equations.

The most commonly used numerical method is the leapfrog integration algorithm
[111] that begins by sampling a value for the momentum variable $\rho$
according to equation 6.11 and then proceeds by applying $L$ times the
following steps,

$\begin{split}\rho&\leftarrow\rho-\frac{\epsilon}{2}\frac{\partial
V}{\partial\theta}\,,\\\ \theta&\leftarrow\theta+\epsilon M^{-1}\rho\,,\\\
\rho&\leftarrow\rho-\frac{\epsilon}{2}\frac{\partial
V}{\partial\theta}\,,\end{split}$ (6.17)

where $\epsilon$ is the integration step size that determines the smallest
time interval. The length of the trajectory will then be $\epsilon L$ and the
new state of the system is denoted as $(\rho^{\prime},\theta^{\prime})$. The
numerical error introduced into the calculation by the leapfrog algorithm is
of the order of $\epsilon^{3}$ per step and $\epsilon^{2}$ globally [111].

##### 28.5 Metropolis acceptance criterion

If the leapfrog algorithm were perfect and did not introduce any numerical
error, we would not have to do anything more than re–sample the momentum
variable every $L$ integration steps. However, the leapfrog integrator is far
from this which means that we need to account for the numerical error that it
introduces before it accumulates. To this end, we only accept and add the new
state $(\rho^{\prime},\theta^{\prime})$ into the Markov chain with probability

$\alpha(\theta^{\prime},\theta)=\min\left(1,\frac{p(\rho^{\prime},\theta^{\prime})}{p(\rho,\theta)}\right)\,,$
(6.18)

and reject it otherwise by adding $(\rho,\theta)$ into the chain. Equation
6.18 is simply the Metropolis acceptance probability for HMC. Therefore, we
see that HMC is essentially a case of Metropolis–Hastings with symmetric
proposal distribution in the augmented state space of $(\rho,\theta)$.

Algorithm 8 Hamiltonian Monte Carlo

0: initial state $\theta_{1}$, potential energy $V(\theta)=-\log p(\theta)$ up
to an additive constant, kinetic energy definition
$T(\rho)=\rho^{T}M^{-1}\rho/2$, number of leapfrog steps $L$, integration step
size $\epsilon$, and number of iterations $N$

0: Markov chain $\theta_{1},\dots,\theta_{N}$ that has $p(\theta)$ as its
equilibrium distribution

1: Draw momentum variable $\rho_{t}\sim\mathcal{N}(0,M)$

2: Set proposed state $\theta^{\prime}\leftarrow\theta_{t}$

3: Set proposed momentum $\rho^{\prime}\leftarrow\rho$

4: for $=1$ to $L$ do

5: Update momentum
$\rho^{\prime}\leftarrow\rho^{\prime}-\frac{\epsilon}{2}\frac{\partial
V}{\partial\theta}$

6: Update position $\theta^{\prime}\leftarrow\theta^{\prime}+\epsilon
M^{-1}\rho^{\prime}$

7: Update momentum
$\rho^{\prime}\leftarrow\rho^{\prime}-\frac{\epsilon}{2}\frac{\partial
V}{\partial\theta}$

8: end for

9: Reverse momentum $\rho^{\prime}\leftarrow-\rho^{\prime}$

10: Compute acceptance probability
$\alpha=\min\left(1,\exp\left[V(\theta)-V(\theta^{\prime})+T(\theta)-T(\theta^{\prime})\right]\right)$

11: Draw uniform number $u\sim\mathcal{U}(0,1)$

12: if $u<\alpha$ then

13: Accept proposed state and set $\theta_{t+1}\leftarrow\theta^{\prime}$

14: else

15: Reject proposed state and set $\theta_{t+1}\leftarrow\theta$

16: end if

##### 28.6 Performance and tuning

The sampling performance of HMC is very sensitive to its tuning [135, 84] and
many efforts have been made to develop heuristics and automated tuning
procedures for the two hyperparameters, $\epsilon$ and $L$, that the method
relies upon. The step size $\epsilon$ can be adaptively tuned by trying to
match the observed acceptance rate to the theoretically optimal value of
$0.65$. Tuning the number of steps $L$ is more cumbersome in practice. In
principle, $L$ can be tuned by minimising the autocorrelation time of the
Markov chain. In practice this requires running multiple preliminary runs with
different values of $L$ in order to determine the most efficient one. For this
reason, other approaches, such as Empirical HMC [187] and the No U-Turn
Sampler (NUTS) [84], have been proposed that automate the use of HMC for many
applications.

### 7 Ensemble MCMC methods

As for me, I am tormented with an everlasting itch for things remote.
I love to sail forbidden seas, and land on barbarous coasts.
— Herman Melville, Moby–Dick or, the Whale

In order to avoid issues caused by multimodality or the need for tuning the
proposal distribution, ensemble MCMC methods rely on an ensemble of parallel
samplers, often called walkers, that sample from an extended probability
distribution. A common way to construct such an extended probability
distribution is using the product density,

$\pi(\\{\theta_{k}\\}_{k=1}^{K})=\prod_{k=1}^{K}p_{k}(\theta_{k})\,,$ (7.1)

where $p_{k}(\theta_{k})$ are the individual densities, one of which can
correspond to the target distribution of interest (e.g. the posterior), and
$K$ is the number of walkers. It is important to note here that $\theta_{k}$
is not the $k$–th component of a vector, but a $D$–dimensional vector itself.

The simplest product density that we can construct based on equation 7.1 is to
assume that $p_{k}(\theta_{k})=p(\theta_{k})$ for all $k$, meaning that the
product density is just the product of $K$ identical copies of the target
distribution $p(\theta_{k})$. A natural question to ask is then why would
anyone want to do this? Why sample $K$ copies of the same distribution instead
of just one? The answer is that the walkers sampling each copy do not have to
be independent of each other and instead are allowed to exchange information
about their current state. For instance, the proposal distribution for a
single walker can depend on the current positions of the rest of the walkers
in the ensemble. This allows for effective proposals that take into account
the relative length–scales and positions of the modes of the target
distribution.

Of course, other product densities, that do not rely on the simplifying
assumption that $p_{k}(\theta_{k})=p(\theta_{k})$ for all $k$, can also be
defined as we will see in the case of the parallel tempering algorithm in
Section 32. In those cases, the goal is not usually to construct effective
proposal distribution but rather to deal with the challenge of multimodality.

#### 29 Gaussian ensemble

Perhaps the simplest way to construct an ensemble MCMC method that limits the
requirement for tuning, to some extent, its proposal distribution is the
Gaussian ensemble (GE) algorithm. GE uses an ensemble of $K$ walkers that
target a product density of the form of equation 7.1, where all copies
$p_{k}(\theta_{k})$ are identical and correspond to the target distribution of
interest (e.g. posterior), and the proposal distribution of each walker is
simply a normal distribution informed by the positions of the rest of the
walkers in the ensemble [165].

Figure 7.1: Illustration of the Gaussian ensemble MCMC method. A new state
$\theta_{k}^{\prime}$ is proposed in the vicinity of the position $\theta_{k}$
of the walker that is updated using an rescaled version of the sample
covariance matrix of the rest of the walkers (i.e. excluding $\theta_{k}$) for
the normal proposal distribution.

In particular, in a given iteration $t$ of the method, the algorithm performs
a loop over the $K$ walkers updating each walker in turn. A new position
$\tilde{\theta}_{k}$ is proposed from a normal distribution,

$\theta_{k}^{\prime}\sim\mathcal{N}\left(\theta|\theta_{k}^{(t-1)},\gamma\Sigma_{-k}\right)\,,$
(7.2)

centred on the current state $\theta_{k}^{(t-1)}$ of the $k$–th walker and
$\gamma$ is a multiplying factor used to scale the covariance matrix in order
to achieve the optimal acceptance rate (e.g. $\gamma=2.38^{2}/D$). The
covariance matrix $\Sigma_{-k}$ of the proposal distribution is simply the
sample covariance estimated using the positions of the ensemble
$\\{\theta_{1}^{(t)},\allowbreak\dots,\allowbreak\theta_{k-1}^{(t)},\allowbreak\theta_{k+1}^{(t-1)},\allowbreak\dots,\allowbreak\theta_{K}^{(t-1)}\\}$
which excludes the $k$–th walker. It is important to notice also that all the
walkers up to and excluding the $k$–th have already been updated and it is
their updated positions that are used to compute the proposal covariance. This
is essentially a Metropolis–within–Gibbs scheme in disguise. The new point
$\tilde{\theta}_{k}$ is then accepted or rejected based on the usual
Metropolis criterion and the process continuous with the next walker until all
of them have been updated.

Algorithm 9 Gaussian ensemble

0: initial state for the ensemble
$\theta^{(1)}=(\theta_{1}^{(1)},\dots,\theta_{K}^{(1)})$, (unnormalised)
target density $f(\theta)\propto p(\theta)$, covariance scaling factor (e.g.
$\gamma=2.38^{2}/D$), and number of iterations $N$

0: Markov chain $\theta_{1},\dots,\theta_{N}$ that has $p(\theta)$ as its
equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $K$ do

3: Compute ensemble mean $\mu_{-k}^{(t)}=\frac{1}{K-1}\sum_{i\neq
k}\theta_{i}^{(t)}$ excluding the $k$th state

4: Compute ensemble covariance matrix
$\Sigma_{-k}^{(t)}=\frac{1}{K-1}\sum_{i\neq
k}(\theta_{i}^{(t)}-\mu_{-k}^{(t)})(\theta_{i}^{(t)}-\mu_{-k}^{(t)})^{T}$
excluding the $k$th state

5: Draw proposal
$\theta_{k}^{\prime}\sim\mathcal{N}(\theta_{k}^{(t)},\gamma\Sigma_{-k}^{(t)})$

6: Compute acceptance probability
$\alpha_{k}=\min\left(1,f(\theta_{k}^{\prime})/f(\theta_{k}^{(t)})\right)$

7: Draw uniform number $u\sim\mathcal{U}(0,1)$

8: if $u<\alpha_{k}$ then

9: Accept proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{\prime}$

10: else

11: Reject proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{(t)}$

12: end if

13: end for

14: end for

GE solves the problem of tuning the proposal, up to the scaling factor
$\gamma$ of the covariance matrix, but still assumes a Gaussian proposal. This
means that we do not expect that GE will perform better than $K$ parallel
well–tuned Random–walk Metropolis samplers. As we will discuss in the next
couple of sections, there are ways to relax this limitation and allow for more
flexible proposals. Last but not least, the estimation of the proposal
covariance matrix requires that the absolute minimum size of the ensemble to
be $D+1$ for the covariance to be non–singular.

#### 30 Affine–invariant stretch move

The affine–invariant ensemble sampler and in particular the stretch move
introduced by [72] is perhaps the most popular ensemble MCMC method in the
astronomical literature, made available in the Python implementation emcee
[50]. The stretch move algorithm relaxes the limitation of the Gaussian
proposal and instead updates each walker in turn along the direction of a
different uniformly selected walker sampled from the rest of the ensemble. As
we will discuss this change introduces both benefits and challenges.

Figure 7.2: Illustration of the affine–invariant stretch move. The selected
walker $\theta_{k}$ is moved to its new position $\theta_{k}^{\prime}$ along
the line defined by $\theta_{j}$ and $\theta_{k}$. $\theta_{j}$ is a walker
that is uniformly selected from the rest of the ensemble (i.e. excluding
$\theta_{k}$).

In particular, in a given iteration $t$ of the method, the algorithm performs
a loop over the $K$ walkers updating each one in turn. In the so–called
stretch move, we move a walker $\theta_{k}$ using a uniformly selected walker
$\theta_{j}$ from the complementary ensemble
$S_{-k}=\\{\theta_{1}^{(t)},\allowbreak\dots,\allowbreak\theta_{k-1}^{(t)},\allowbreak\theta_{k+1}^{(t-1)},\allowbreak\dots,\allowbreak\theta_{K}^{(t-1)}\\}$
that excludes $\theta_{k}$. The $\theta_{j}$ walker acts as an anchor point
for the move that consists of a proposal of the form

$\theta_{k}^{\prime}=\zeta\theta_{k}+(1-\zeta)\theta_{j}\,,$ (7.3)

where $\zeta$ is a scaling variable with a probability density $g$ that
satisfies the symmetry condition,

$g\left(\frac{1}{\zeta}\right)=\zeta g(\zeta)\,,$ (7.4)

such that the move expressed by equation 7.3 is symmetric in the Metropolis
sense. A particular density that obeys this condition is

$g(\zeta)\propto\begin{cases}\frac{1}{\sqrt{\zeta}}&\mathrm{if}\;\zeta\in\left[\frac{1}{\alpha},\alpha\right]\,,\\\
0&\mathrm{otherwise}\,,\end{cases}$ (7.5)

where $\alpha>1$ is a parameter that can be tuned to enhance the performance.
The default value is usually set to $\alpha=2$. The new state
$\theta_{k}^{\prime}$ is then accepted with Metropolis probability,

$\alpha(\theta_{k}^{\prime},\theta_{k})=\min\left(1,\zeta^{D-1}\frac{p(\theta_{k}^{\prime})}{p(\theta_{k})}\right)\,,$
(7.6)

where the $\zeta^{D-1}$ comes from the fact that the update takes place along
a straight line. The process is then repeated for the next walker, until all
the walkers are updated for the current iteration $t$ before the algorithm
moves to its next iteration.

Algorithm 10 Affine–invariant stretch move

0: initial state for the ensemble
$\theta^{(1)}=(\theta_{1}^{(1)},\dots,\theta_{K}^{(1)})$, (unnormalised)
target density $f(\theta)\propto p(\theta)$, and number of iterations $N$

0: Markov chain $\theta_{1},\dots,\theta_{N}$ that has $p(\theta)$ as its
equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $K$ do

3: Draw a walker $\theta_{j}$ from the complementary ensemble
$S_{-k}=\left\\{\theta_{1}^{(t+1)},\dots,\theta_{k-1}^{(t+1)},\theta_{k+1}^{(t)},\dots,\theta_{K}^{(t)}\right\\}$

4: Draw random number $\zeta\sim g(\zeta)$

5: Compute proposed state
$\theta_{k}^{\prime}\leftarrow\zeta\theta_{k}+(1-\zeta)\theta_{j}$

6: Compute acceptance probability
$\alpha_{k}=\min\left(1,\zeta^{D-1}f(\theta_{k}^{\prime})/f(\theta_{k}^{(t)})\right)$

7: Draw uniform number $u\sim\mathcal{U}(0,1)$

8: if $u<\alpha_{k}$ then

9: Accept proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{\prime}$

10: else

11: Reject proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{(t)}$

12: end if

13: end for

14: end for

One of the strict requirements of this method is the minimum number of walkers
to be $D+1$ for it to be ergodic and avoid the risk of walkers getting trapped
in some hyper–plane of lower than $D$ dimensions. Practically, the actual
number of walkers required is much larger as it determines the plethora of
possible new directions along which updates take place in each iteration. In
this sense, the initial positions of the walkers and the number of them are
the only free hyperparameters of this method. A great benefit of this method
is that it is affine–invariant, that is, its performance is insensitive to any
linear correlations between the parameters of the target distribution. As the
astronomical community has witnessed during the past few years, this offers a
great advantage over other methods.

#### 31 Differential evolution

Another ensemble method in the spirit of the stretch move is the differential
evolution MCMC [169, 170]. Unlike the stretch move that requires another
single walker to act as an anchor point for a proposal, differential evolution
involves two. We will discuss shortly how this difference can affect the
performance and alter the characteristics of the method.

Figure 7.3: Illustration of the differential evolution Monte Carlo. The
selected walker $\theta_{k}$ is moved to its new position
$\theta_{k}^{\prime}$ parallel to the line defined by $\theta_{i}$ and
$\theta_{j}$. The latter are two walkers that are uniformly selected from the
rest of the ensemble (i.e. excluding $\theta_{k}$).

An update of the ensemble works as follows: the algorithm performs a loop over
the $K$ walkers updating each one in turn. Assuming that the current walker to
be updated is $\theta_{k}$, the algorithm uniformly selects two walkers
(without replacement), $\theta_{i}$ and $\theta_{j}$, from the complementary
ensemble
$S_{-k}=\\{\theta_{1}^{(t)},\allowbreak\dots,\allowbreak\theta_{k-1}^{(t)},\allowbreak\theta_{k+1}^{(t-1)},\allowbreak\dots,\allowbreak\theta_{K}^{(t-1)}\\}$
that excludes $\theta_{k}$. The vector $\theta_{i}-\theta_{j}$ connecting the
two auxiliary walkers defines the direction along which a move is proposed.
The move consists of a proposal of the form

$\theta_{k}^{\prime}=\theta_{k}+\gamma\times\left(\theta_{i}-\theta_{j}\right)+\epsilon\,,$
(7.7)

where $\gamma$ is a non–zero scaling factor and
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$ is some optional Gaussian noise. The
value of $\gamma$ determines the scale of the proposal. Its default value is
often set to $\gamma=2.38/\sqrt{2D}$ which results in the optimal acceptance
rate (i.e. $23.4\%$) for normal target distributions. In practice, one can
adapt $\gamma$ using some diminishing adaptation scheme during the run. The
proposed update of equation 7.7 is then accepted with Metropolis acceptance
probability

$\alpha(\theta^{\prime},\theta)=\min\left(1,\frac{p(\theta^{\prime})}{p(\theta)}\right)\,.$
(7.8)

Algorithm 11 Differential evolution

0: initial state for the ensemble
$\theta^{(1)}=(\theta_{1}^{(1)},\dots,\theta_{K}^{(1)})$, (unnormalised)
target density $f(\theta)\propto p(\theta)$, proposal scale parameter (e.g.
$\gamma=2.38/\sqrt{2D}$), optional Gaussian noise standard deviation (e.g.
$\sigma=10^{-3})$, and number of iterations $N$

0: Markov chain $\theta_{1},\dots,\theta_{N}$ that has $p(\theta)$ as its
equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $K$ do

3: Draw walkers $\theta_{i}$ and $\theta_{j}$ without replacement from the
complementary ensemble
$S_{-k}=\left\\{\theta_{1}^{(t+1)},\dots,\theta_{k-1}^{(t+1)},\theta_{k+1}^{(t)},\dots,\theta_{K}^{(t)}\right\\}$

4: Draw random noise $\epsilon\sim\mathcal{N}(0,\sigma^{2})$

5: Compute proposed state
$\theta_{k}^{\prime}\leftarrow\theta_{k}+\gamma(\theta_{i}-\theta_{j})+\epsilon$

6: Compute acceptance probability
$\alpha_{k}=\min\left(1,f(\theta_{k}^{\prime})/f(\theta_{k}^{(t)})\right)$

7: Draw uniform number $u\sim\mathcal{U}(0,1)$

8: if $u<\alpha_{k}$ then

9: Accept proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{\prime}$

10: else

11: Reject proposed state and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{(t)}$

12: end if

13: end for

14: end for

The advantage of differential evolution over the stretch moves comes down to
the flexibility of their proposals. The direction along which a walker moves
in the context of the stretch moves is determined by a single walker. This
means that at any given iteration, the number of equally possible directions
is $K-1$. On the other hand, differential evolution moves each walker along a
direction defined by two walkers. This implies that the total number of
possible directions is given by the binomial combination $\binom{K-1}{2}$. The
latter increases much faster with the number size of the ensemble $K$ than the
former, offering a larger variety of possible trajectories for the walkers. In
other words, differential evolution is expected to perform better even with a
lower number of walkers.

#### 32 Parallel tempering

So far we have only discussed ensemble methods that target a trivial product
density given by the product of $K$ copies of the target distribution as shown
in equation 7.1. The main rationale for attempting to do this was to reduce
the tuning requirements of MCMC. If we focus on addressing the challenge of
multimodality, that is, the existence of multiple peaks in the target
distribution, then we have to introduce a different product density as the
extended target distribution.

One such choice is,

$p^{*}\left(\\{\theta_{k}\\}_{k=1}^{K}\right)=\prod_{k=1}^{K}p_{k}(\theta_{k})\,,$
(7.9)

where

$p_{k}(\theta_{k})\propto
p^{\beta_{k}}(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\,,$ (7.10)

is the annealed or tempered posterior that offers a simple interpolation
between the prior $p(\theta|\mathcal{M})$ and the unnormalised posterior
density $p(d|\theta,\mathcal{M})\allowbreak p(\theta|\mathcal{M})$ for
different monotonically–increasing values of $\beta_{k}\in[0,1]$. In the limit
that $\beta_{k}=1$ for all values of $k$, equation 7.9 reduces to the usual
product density of equation 7.1.

Figure 7.4: Illustration of the gradual tempering performed in the posterior
distribution. The prior distribution corresponds to $\beta\rightarrow 0$ and
the posterior is recovered as $\beta\rightarrow 1$.

The method of parallel tempering (PT) [40], also known as replica exchange
Monte Carlo (REMC) [167, 87] or Metropolis–coupled Markov chain Monte Carlo
(MC3) [66], relies on $K$ parallel Markov chains, each one targeting a
different tempered density. The $\beta_{k}$ values are usually chosen a priori
using a heuristic rule (e.g. $\beta_{k}=(k-1)^{3}/(K-1)^{3}$), or are set
adaptively during the run using some diminishing adaptation scheme. The choice
of MCMC method used for each different $\beta_{k}$ is completely arbitrary and
it can be anything from simple Random–walk Metropolis to Hamiltonian Monte
Carlo or even an ensemble MCMC method.

So far, PT might look very similar to a parallel version of the sequential
simulated annealing method in which $\beta_{k}=1/T_{k}$ has the role of the
inverse temperature. The crucial difference that makes PT so powerful is the
fact that one can also perform between–chain exchange moves. Either
periodically (e.g. once every $10$ steps) or randomly (e.g. with probability
$10\%$) a swap can take place between two states $\theta_{i}$ and $\theta_{j}$
that belong to different tempered posteriors (i.e. $\beta_{i}\neq\beta_{j}$).
The reason that exchange/swap moves are desirable is that they enable the
transfer of information from states of low $\beta$ to those of higher $\beta$.

Figure 7.5: Illustration of the parallel tempering swaps performed between
adjacent temperature levels.

To understand how to perform a swap in practice let us consider the extended
state,

$\\{\theta_{k}\\}_{k=1}^{K}=\\{\theta_{1},\dots,\theta_{i},\dots,\theta_{j},\dots,\theta_{K}\\}\,,$
(7.11)

prior to the swap, where $\theta_{i}$ and $\theta_{j}$ are the two states that
we want to exchange. This means that the proposed new state will be,

$\\{\tilde{\theta}_{k}\\}_{k=1}^{K}=\\{\theta_{1},\dots,\theta_{j},\dots,\theta_{i},\dots,\theta_{K}\\}\,.$
(7.12)

Notice that the rest of the states, with the exception of $\theta_{i}$ and
$\theta_{j}$, are left unaffected by this exchange proposal. The Metropolis
acceptance probability for this proposal is,

$\begin{split}\alpha_{ij}&=\min\left(1,\frac{p^{*}\left(\\{\tilde{\theta}_{k}\\}_{k=1}^{K}\right)}{p^{*}\left(\\{\theta_{k}\\}_{k=1}^{K}\right)}\right)\\\
&=\min\left(1,\frac{p_{i}(\theta_{j})p_{j}(\theta_{i})}{p_{i}(\theta_{i})p_{j}(\theta_{j})}\right)\\\
&=\min\left(1,\frac{p^{\beta_{i}}(d|\theta_{j},\mathcal{M})p(\theta_{j}|\mathcal{M})p^{\beta_{j}}(d|\theta_{i},\mathcal{M})p(\theta_{i}|\mathcal{M})}{p^{\beta_{i}}(d|\theta_{i},\mathcal{M})p(\theta_{i}|\mathcal{M})p^{\beta_{j}}(d|\theta_{j},\mathcal{M})p(\theta_{j}|\mathcal{M})}\right)\\\
&=\min\left(1,\frac{p^{\beta_{i}}(d|\theta_{j},\mathcal{M})p^{\beta_{j}}(d|\theta_{i},\mathcal{M})}{p^{\beta_{i}}(d|\theta_{i},\mathcal{M})p^{\beta_{j}}(d|\theta_{j},\mathcal{M})}\right)\\\
&=\min\left[1,\left(\frac{p(d|\theta_{i},\mathcal{M})}{p(d|\theta_{j},\mathcal{M})}\right)^{(\beta_{j}-\beta_{i})}\right]\,,\end{split}$
(7.13)

The chains are usually chosen to be in adjacent $\beta_{k}$ levels (i.e.
$i=j-1$) such that the overlap between the typical sets of $p_{i}(\theta_{i})$
and $p_{j}(\theta_{j})$ is large, leading to high acceptance probabilities.
[6] estimated that the optimal acceptance rate is $23.4\%$. The exchange
updates are typically performed after the local MCMC updates are completed in
all $\beta$ levels for a given iteration. Furthermore, there are different
strategies for proposing swaps between adjacent temperature levels [114]. One
option is to randomly select a pair of adjacent temperature levels per
iteration. Another strategy involves proposing to swap all adjacent pairs
starting from the lowest or highest $\beta$ level and progressively moving
towards the other end of the ladder. Finally, strategies that involve two
steps, for instance, proposing to swap all even pairs in one iteration and all
odd pairs in the next, have also been suggested in the literature [114].

Algorithm 12 Parallel tempering

0: initial state for the ensemble
$\theta^{(1)}=(\theta_{1}^{(1)},\dots,\theta_{K}^{(1)})$, prior probability
density $\pi(\theta)\equiv p(\theta|\mathcal{M})$, likelihood function
$\mathcal{L}(\theta)\equiv p(d|\theta,\mathcal{M})$, temperature ladder (e.g.
$\beta_{k}=(k-1)^{3}/(K-1)^{3}$), local MCMC kernel
$\theta^{\prime}\leftarrow\mathcal{K}(\theta;f(\theta))$ (e.g. a single
random–walk Metropolis update), and number of iterations $N$

0: $K$ Markov chains that each has
$p_{t}(\theta)\propto\pi(\theta)\mathcal{L}(\theta)^{\beta_{k}}$ as its
equilibrium distribution

1: for $t=1$ to $N$ do

2: for $k=1$ to $K$ do

3: Update state using local MCMC update
$\theta_{k}^{\prime}\leftarrow\mathcal{K}(\theta_{k}^{(t)};\pi(\theta)\mathcal{L}(\theta)^{\beta_{k}})$

4: end for

5: Draw random value of $k$ uniformly $k\sim\mathcal{U}(1,K-1)$

6: Compute acceptance probability
$\alpha_{k}=\min\left(1,\left[\frac{\mathcal{L}\left(\theta_{k}^{(t+1)}\right)}{\mathcal{L}\left(\theta_{k+1}^{(t+1)}\right)}\right]^{\beta_{k+1}-\beta_{k}}\right)$

7: Draw uniform number $u\sim\mathcal{U}(0,1)$

8: if $u<\alpha_{k}$ then

9: Accept proposed swap and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k+1}^{\prime}$ and
$\theta_{k+1}^{(t+1)}\leftarrow\theta_{k}^{\prime}$

10: else

11: Reject proposed swap and set
$\theta_{k}^{(t+1)}\leftarrow\theta_{k}^{\prime}$ and
$\theta_{k+1}^{(t+1)}\leftarrow\theta_{k+1}^{\prime}$

12: end if

13: end for

### 8 Evidence and Bayes factor computation

There is nothing more deceptive than an obvious fact.
— Arthur Conan Doyle

#### 33 Naive Monte Carlo estimator

The simplest estimator for the evidence we can construct is just the
expectation value of the likelihood function with respect to the prior
distribution [78, 145]. The, so–called, Naive Monte Carlo (NMC) estimator can
be computed as the sum

$\hat{p}_{\mathrm{NMC}}(d|\mathcal{M})=\frac{1}{n}\sum_{i=1}^{n}p(d|\theta_{i},\mathcal{M})\,,\quad\mathrm{with}\quad\theta_{i}\sim
p(\theta|\mathcal{M})\,.$ (8.1)

Although simple and unbiased, this approach can become extremely inefficient
and result in a high variance in higher dimensions as the probability mass
concentrates in the typical set that occupies a negligible fraction of the
prior volume [137]. For this reason, this technique is only recommended for
low–dimensional problems (i.e. $D\leq 3$).

#### 34 Importance sampling estimator

A more general strategy for the unbiased estimation of the evidence is
importance sampling using samples from an auxiliary distribution $q(\theta)$.
A simple estimator can then be constructed as,

$\hat{p}_{\mathrm{IS}}(d|\mathcal{M})=\frac{1}{n}\sum_{i=1}^{n}\frac{p(\theta_{i}|d,\mathcal{M})}{q(\theta_{i})}\,,\quad\mathrm{with}\quad\theta_{i}\sim
q(\theta)\,.$ (8.2)

These estimators share the same difficulty as most methods based upon
importance sampling, that is, a large overlap between the typical set of the
proposal and posterior distribution must be achieved for the method to be
effective. Constructing effective proposal distributions becomes increasingly
unmanageable as the number of dimensions increases and thus the application of
this method on its own is limited to low dimensions. Finally, the importance
sampling estimator reduces to the NMC one when the proposal distribution is
chosen to be the prior.

#### 35 Harmonic mean estimator

The harmonic mean (HM) estimator is another variation of the importance
sampling estimator in which the posterior is used as the proposal and the
prior as the target distribution [137]. This suggests the following estimator,

$\hat{p}_{\mathrm{HM}}(d|\mathcal{M})=\frac{1}{\frac{1}{n}\sum_{i=1}^{n}\frac{1}{p(d|\theta_{i}\mathcal{M})}}\,,\quad\mathrm{with}\quad\theta_{i}\sim
p(\theta|d,\mathcal{M})\,.$ (8.3)

The possible occurrence of samples with small likelihood value renders the
variance of this estimator infinite [134]. This pathology can be addressed by
using a mixture $q(\theta)=\delta
p(\theta|\mathcal{M})+(1-\delta)p(\theta|d,\mathcal{M})$ between the prior and
the posterior as the proposal distribution, where $\delta$ is very small (e.g.
$\delta=0.05$). The resulting method is then called the stabilised harmonic
mean (SHM) estimator [137].

#### 36 Laplace estimator

As discussed in detail in Section 20, for a sufficiently Gaussian target
distribution $p(\theta)$ we can use the Laplace approximation, that is, a
second order expansion around the mode, to estimate expectation values [173].
Assuming that the target distribution is the unnormalised posterior
$p(d|\theta\mathcal{M})p(\theta|\mathcal{M})$, the Gaussian approximation’s
mean is given by,

$\mu=\underset{\theta}{\mathrm{arg\,max}}\leavevmode\nobreak\
\left[p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\right]\,,$ (8.4)

following equation 4.23, and the inverse covariance is given by,

$\left(\Sigma^{-1}\right)_{ij}=-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\left[\log
p(d|\theta,\mathcal{M})+\log p(\theta|\mathcal{M})\right]\,,$ (8.5)

following equation 4.24. Then, the model evidence is approximated by the
normalising constant of the Gaussian, or in other words,

$\begin{split}\hat{p}_{L}(d|\mathcal{M})&=\int
e^{-\frac{1}{2}(\theta-\mu)^{T}\Sigma^{-1}(\theta-\mu)}d\theta\\\
&=(2\pi)^{D/2}\det(\Sigma)^{1/2}p(d|\theta=\mu,\mathcal{M})p(\theta=\mu|\mathcal{M})\,.\end{split}$
(8.6)

As with any method, this result is only as good as the assumptions entering
its calculation. The closer the posterior resembles a normal distribution, the
better the outcome of the Laplace estimator will be.

#### 37 Bridge sampling

Originally, [120] introduced bridge sampling (BS) as a way to directly
estimate the Bayes factor of two models, $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$. However, in this section we present a version of BS that
targets the model evidence of a single model $\mathcal{M}$. BS follows from
the basic identity,

$1=\frac{\int
p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\alpha(\theta)q(\theta)}{\int
p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\alpha(\theta)q(\theta)}\,,$ (8.7)

where $q(\theta)$ is the proposal distribution and $\alpha(\theta)$ is the
so–called bridge function the support of which encompasses that of both the
target posterior and of the proposal distribution.

Multiplying both sides of equation 8.7 with the model evidence
$p(d|\mathcal{M})$ results in

$\begin{split}p(d|\mathcal{M})&=\frac{\int
p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\alpha(\theta)q(\theta)}{\int\frac{p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})}{p(d|\mathcal{M})}\alpha(\theta)q(\theta)}\\\
&=\frac{\int
p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\alpha(\theta)q(\theta)}{\int
p(\theta|d,\mathcal{M})\alpha(\theta)q(\theta)}\,,\end{split}$ (8.8)

which can be written as,

$p(d|\mathcal{M})=\frac{\mathbb{E}_{q(\theta)}\left[p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\alpha(\theta)\right]}{\mathbb{E}_{p(\theta|d,\mathcal{M})}\left[\alpha(\theta)q(\theta)\right]}\,,$
(8.9)

in terms of expectation values. The model evidence can then be approximated
as,

$\hat{p}(d|\mathcal{M})=\frac{\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}p(d|\tilde{\theta}_{i},\mathcal{M})p(\tilde{\theta}_{i}|\mathcal{M})\alpha(\tilde{\theta}_{i})}{\frac{1}{n_{1}}\sum_{j=1}^{n_{1}}\alpha(\theta_{j}^{*})q(\theta_{j}^{*})}\,,$
(8.10)

where $\tilde{\theta}_{i}$ are samples from the proposal distribution,

$\tilde{\theta}_{i}\sim q(\theta)\,,$ (8.11)

and $\theta_{j}^{*}$ are samples from the posterior distribution,

$\theta_{j}^{*}\sim p(\theta|d,\mathcal{M})\,.$ (8.12)

It is clear from the above discussion that BS relies on samples from both the
proposal distribution $q(\theta)$, which plays the role of an importance
density, and the posterior distribution $p(\theta|d,\mathcal{M})$. Often, the
proposal distribution is some distribution that is easy to sample from and its
typical set has a large overlap with the one of the posterior distribution. A
common proposal used in practice is a normal distribution with its first two
moments matching those of the posterior distribution.

Although highly arbitrary, the choice of the bridge function $\alpha(\theta)$
can have a significant impact on the precision of the method for a given
proposal distribution. For instance, setting $\alpha(\theta)=[q(\theta)]^{-1}$
the BS estimator reduces to the naive Monte Carlo estimator, whereas setting
$\alpha(\theta)=[p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})q(\theta)]^{-1}$
leads to the harmonic mean estimator. [120] showed that the optimal bridge
function, that is, the one that minimises the mean–square–error, is,

$\alpha(\theta)=\frac{C}{s_{1}p(d|\theta,\mathcal{M})p(\theta|\mathcal{M})+s_{2}p(d|\mathcal{M})q(\theta)}\,,$
(8.13)

where $s_{1}=n_{1}/(n_{1}+n_{2})$ and $s_{2}=n_{2}/(n_{1}+n_{2})$ and $C$ is a
constant that cancels out and its value does not affect the outcome in any
way. The bridge function of equation 8.13 depends on the model evidence
$p(d|\mathcal{M})$, the same quantity that we are trying to approximate. We
can resolve this issue by employing an iterative scheme,

$\hat{p}_{\mathrm{BS}}^{(t+1)}(d|\mathcal{M})=\frac{\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\frac{p(d|\tilde{\theta}_{i},\mathcal{M})p(\tilde{\theta}_{i}|\mathcal{M})}{s_{1}p(d|\tilde{\theta}_{i},\mathcal{M})p(\tilde{\theta}_{i}|\mathcal{M})+s_{2}\hat{p}_{\mathrm{BS}}^{(t)}(d|\mathcal{M})q(\tilde{\theta}_{i})}}{\frac{1}{n_{1}}\sum_{j=1}^{n_{1}}\frac{q(\theta_{j}^{*})}{s_{1}p(d|\theta_{j}^{*},\mathcal{M})p(\theta_{j}^{*}|\mathcal{M})+s_{2}\hat{p}_{\mathrm{BS}}^{(t)}(d|\mathcal{M})q(\theta_{j}^{*})}}\,,$
(8.14)

starting from some initial guess of the value of the model evidence
$\hat{p}_{\mathrm{BS}}^{(0)}(d|\mathcal{M})$ and keep updating it until the
estimate has converged for some arbitrary tolerance level. Rearranging the
terms on the right hand side, the aforementioned estimator can be written in
the simpler form

$\hat{p}_{\mathrm{BS}}^{(t+1)}(d|\mathcal{M})=\frac{\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\frac{\ell_{2,i}}{s_{1}\ell_{2,i}+s_{2}\hat{p}_{\mathrm{BS}}^{(t)}(d|\mathcal{M})}}{\frac{1}{n_{1}}\sum_{j=1}^{n_{1}}\frac{1}{s_{1}\ell_{1,j}+s_{2}\hat{p}_{\mathrm{BS}}^{(t)}(d|\mathcal{M})}}\,,$
(8.15)

where we have defined

$\ell_{1,j}=\frac{p(d|\theta_{j}^{*},\mathcal{M})p(\theta_{j}^{*}|\mathcal{M})}{q(\theta_{j}^{*})}\,,$
(8.16)

and

$\ell_{2,i}=\frac{p(d|\tilde{\theta}_{i},\mathcal{M})p(\tilde{\theta}_{i}|\mathcal{M})}{q(\tilde{\theta}_{i})}\,.$
(8.17)

Furthermore, the numerical stability of equation 8.15 can be improved and
overflow issues avoided if we define

$\hat{p}_{\mathrm{BS}}^{(t)}(d|\mathcal{M})=\hat{r}^{(t)}\exp(\ell^{*})\,,$
(8.18)

and use the iterative formula

$\hat{r}^{(t+1)}=\frac{\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\frac{\exp[\log(\ell_{2,i})-\ell^{*}]}{s_{1}\exp[\log(\ell_{2,i})-\ell^{*}]+s_{2}\hat{r}^{(t)}}}{\frac{1}{n_{1}}\sum_{j=1}^{n_{1}}\frac{1}{s_{1}\exp[\log(\ell_{1,j})-\ell^{*}]+s_{2}\hat{r}^{(t)}}},$
(8.19)

where $\ell^{*}$ is a constant that we can choose in order to make the sums
numerically tractable, for instance
$\ell^{*}=\mathrm{median}[\log(\ell_{1,j})]$.

Compared to other methods such as importance sampling or the harmonic mean
estimator, BS estimates are more robust in cases in which the overlap between
the typical sets of the proposal and posterior distribution is far from
perfect.

#### 38 Thermodynamic integration

A large body of work in statistical physics is concerned with methods for the
estimation of normalising constants and partition functions in particular. The
method of thermodynamic integration (TI) was developed for exactly this
purpose [61]. [51] studied the particular case in which the normalising
constant that is estimated using TI is the model evidence. To this end, they
introduced the notion of the power posterior,

$p(\theta|d,\beta,\mathcal{M})\propto
p^{\beta}(d|\theta,\mathcal{M})p(\theta|\mathcal{M})\,,$ (8.20)

in which $\beta$ is an auxiliary variable in the interval $[0,1]$. By
construction, the normalising constant of the power posterior is simply,

$p(d|\beta,\mathcal{M})=\int
p^{\beta}(d|\theta,\mathcal{M})p(\theta|\mathcal{M})d\theta\,,$ (8.21)