where $p(d|\beta=1,\mathcal{M})$ is the model evidence and
$p(d|\beta=0,\mathcal{M})$ is the integral over the prior which is simply
equal to $1$. Furthermore, the logarithm of the model evidence is,

$\begin{split}\log
p(d|\mathcal{M})&=\log\left[\frac{p(d|\beta=1,\mathcal{M})}{p(d|\beta=0,\mathcal{M})}\right]\\\
&=\int_{0}^{1}\mathbb{E}_{p(\theta|d,\beta,\mathcal{M})}[\log
p(d|\theta,\beta,\mathcal{M})]d\beta\,,\end{split}$ (8.22)

that is, the integral over $\beta$ of the expectation value of the likelihood
with respect to the posterior for each value of $\beta$. To prove the above
identity we first need to notice that,

$\begin{split}\frac{d\log
p(d|\beta,\mathcal{M})}{d\beta}&=\frac{1}{p(d|\beta,\mathcal{M})}\times\frac{dp(d|\beta,\mathcal{M})}{d\beta}\\\
&=\frac{1}{p(d|\beta,\mathcal{M})}\int\frac{d}{d\beta}p^{\beta}(d|\theta,\mathcal{M})p(\theta|\mathcal{M})d\theta\\\
&=\int\log
p(d|\theta,\mathcal{M})\frac{p^{\beta}(d|\theta,\mathcal{M})p(\theta|\mathcal{M})}{p(d|\beta,\mathcal{M})}d\theta\\\
&=\int\log p(d|\theta,\mathcal{M})p(\theta|d,\beta,\mathcal{M})d\theta\\\
&=\mathbb{E}_{p(\theta|d,\beta,\mathcal{M})}\left[\log
p(d|\theta,\mathcal{M})\right]\,,\end{split}$ (8.23)

Integrating both sides with respect to $\beta$ leads to equation 8.22 and
completes the proof.

Using equation 8.22 to estimate the model evidence often requires the
discretisation of the integral. A sequence
$0=\beta_{1}<\beta_{2}<\dots<\beta_{m}=1$ must be chosen a priori or based on
diminishing adaptation scheme. The model evidence can then be approximated
using the trapezoidal rule,

$\hat{p}_{\mathrm{TI}}(d|\mathcal{M})=\sum_{i=1}^{m}(\beta_{i+1}-\beta_{i})\frac{E_{i}+E_{i+1}}{2}\,,$
(8.24)

where,

$E_{i}=\mathbb{E}_{p(\theta|d,\beta_{i},\mathcal{M})}\left[\log
p(d|\theta,\mathcal{M})\right]\,,$ (8.25)

is the expected likelihood at $\beta_{i}$.

There are two sources of error in the above approximation. The first one is
the Monte Carlo error that originates from the estimation of equation 8.25
using a finite number of samples. The second type has to do with the choice of
the discretisation of $\beta$. [25] showed that the discretisation error
depends on the Kullback–Leibler (KL) divergence between subsequent densities
$p(\theta|d,\beta_{i},\mathcal{M})$ and $p(\theta|d,\beta_{i+1},\mathcal{M})$.
This means that the optimal discretisation sequence of $\beta$ values is the
one minimising the KL divergence between subsequent power posteriors. Of
course, knowing the optimal scheme is a priori hardly ever possible and thus
we must rely on ad hoc choices (e.g. $\beta_{i}=(i-1)^{3}/(m-1)^{3}$) or
diminishing adaptation strategies.

Thermodynamic integration can be combined with many different MCMC methods in
order to estimate the model evidence. Perhaps the simplest one is to run $m$
independent chains, either in parallel or serially, and then estimate the
evidence using equation 8.24 where the expected likelihood of each discrete
$\beta$ value is computed with equation 8.25 for each chain. Furthermore, the
chains do not even have to be independent for this method to work. Lastly, a
parallel tempering approach can be followed as it is often done in practice.

#### 39 Annealed importance sampling

Annealed importance sampling (AIS) is another method that relies on a sequence
of annealed or tempered distributions in order to construct an importance
sampling estimator for the model evidence [132], similarly to simulated
annealing and parallel tempering.

The basic idea is to use MCMC transitions in order to push a collection of $n$
particles through a series of $m$ intermediate distributions

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

where $0=\beta_{1}<\beta_{2}<\dots<\beta_{m}=1$, connecting the prior for
$\beta_{1}=0$ to the posterior for $\beta_{m}=1$. The particles are
initialised by drawing samples from the prior

$\theta_{i}^{(1)}\sim p(\theta|\mathcal{M})\,,$ (8.27)

and assigned (unnormalised) importance weights

$w_{i}^{(1)}=1\,,$ (8.28)

for $i\in\\{1,\dots,n\\}$ the particle index.

A number of $N$ MCMC steps is then performed for each particle before the
value of $\beta$ is updated to the next value in the predefined sequence. The
number $N$ of MCMC steps is chosen such that the Markov chains defined by the
particle trajectories have enough time to reach the stationary distribution.
The critical difference between AIS and simulated annealing is that the
associated importance weights $w_{i}$ are updated during the run every time we
move from one intermediate distribution to the next,

$w_{i}^{(t+1)}=w_{i}^{(t)}\times\frac{p(\theta_{i}|d,\beta_{t+1},\mathcal{M})}{p(\theta_{i}|d,\beta_{t},\mathcal{M})}=w_{i}^{(t)}\times\frac{p^{\beta_{t+1}}(d|\theta_{i},\mathcal{M})}{p^{\beta_{t}}(d|\theta_{i},\mathcal{M})}\,.$
(8.29)

In practice, the logarithm of the weights is used in order to avoid numerical
issues. Once the final distribution (i.e. the posterior) is reached and the
particle weights are updated accordingly, the model evidence can be estimated
as

$\hat{p}_{\mathrm{AIS}}(d|\mathcal{M})=\frac{1}{n}\sum_{i=1}^{n}w_{i}^{(m)}\,.$
(8.30)

Furthermore, the samples $\theta_{i}^{(m)}$ combined with their respective
weights $w_{i}^{(m)}$ can be used to compute arbitrary expectation values

$\mathbb{E}_{p(\theta|d,\mathcal{M})}\left[f(\theta)\right]=\frac{\sum_{i=1}^{n}w_{i}^{(m)}f\left(\theta_{i}^{(m)}\right)}{\sum_{i=1}^{n}w_{i}^{(m)}}\,.$
(8.31)

Assuming that the annealing process is slow enough (i.e. large number $m$ of
$\beta$ levels and number $N$ of MCMC steps) and a large enough collection of
particles is used, then AIS yields unbiased estimates of the model evidence
and weighted posterior samples, even in high dimensions. In the limit that the
number of MCMC steps $N$ goes to zero, the AIS estimator reduces to the usual
importance sampling estimator.

#### 40 Savage–Dickey density ratio

Suppose now that we have two models or hypotheses and their respective
parameters, $\mathcal{M}_{0}:\theta$ and $\mathcal{M}_{1}:\theta,\phi$, such
that $\mathcal{M}_{0}$ is nested inside $\mathcal{M}_{1}$. This means that the
more complex model, $\mathcal{M}_{1}$, is reduced to the simpler one,
$\mathcal{M}_{0}$, for some specific choice of one or more of its parameters,
$\phi=\phi_{0}$. This specific parameter choice is often called a point–null
hypothesis as it is associated with zero probability mass in the context of
the $\mathcal{M}_{1}$ model.

Evaluating the plausibility of this hypothesis can be done by computing the
Bayes factor between the two models. The Savage–Dickey density ratio (SDDR) is
a method that aims to do exactly this that was introduced by [37], [35], [76],
and [36] who in turn attributed the origin of the method to Leonard Jimmie
Savage.

Although the SDDR can be only applied to nested models, it has the advantage
that it is simple to compute, given some posterior samples, without making any
assumptions about the Gaussianity of the posterior distribution. In
particular, the Bayes factor of $\mathcal{M}_{0}$ over $\mathcal{M}_{1}$ is
simply,

$\mathrm{BF}_{01}=\frac{p(d|\mathcal{M}_{0})}{p(d|\mathcal{M}_{1})}=\frac{p(\phi=\phi_{0}|d,\mathcal{M}_{1})}{p(\phi=\phi_{0}|\mathcal{M}_{1})}\,,$
(8.32)

where the numerator of the right–hand–side ratio is just the marginal
posterior of $\phi$ for $\mathcal{M}_{1}$ evaluated at $\phi=\phi_{0}$, and
the denominator is the prior of $\phi$ for $\mathcal{M}_{1}$ evaluated at
$\phi=\phi_{0}$. In other words, the Bayes factor is simply the marginal
posterior to prior ratio for $\mathcal{M}_{1}$ evaluated at $\phi=\phi_{0}$.
This means that only the parameters $\phi$ determine the value of the Bayes
factor, and the nuisance parameters $\theta$, that are common among the two
models, are irrelevant. A schematic representation of SDDR is depicted in
Figure 8.1.

Figure 8.1: The Savage–Dickey density ratio expresses the Bayes factor
$\mathrm{BF}_{01}$ as the ratio of marginal posterior to the prior density at
the point $\phi_{0}$ in which model $M_{1}$ reduces to $M_{0}$.

The validity of this method relies on two conditions. First, that the
likelihood function of $\mathcal{M}_{1}$ has to reduce to that of
$\mathcal{M}_{0}$

$p(d|\theta,\phi=\phi_{0},\mathcal{M}_{1})=p(d|\theta,\mathcal{M}_{0})\,,$
(8.33)

and the same must be true for the prior

$p(\theta|\phi=\phi_{0},\mathcal{M}_{1})=p(\theta|\mathcal{M}_{0})\,.$ (8.34)

The condition of equation 8.34 is satisfied by separable priors,

$p(\theta,\phi|\mathcal{M}_{1})=p(\theta|\mathcal{M}_{1})p(\phi|\mathcal{M}_{1})\,.$
(8.35)

The proof of equation 8.32 is straightforward, starting with,

$\begin{split}p(d|\mathcal{M}_{0})&=\int
p(d|\theta,\mathcal{M}_{0})p(\theta|\mathcal{M}_{0})d\theta\\\ &=\int
p(d|\theta,\phi=\phi_{0},\mathcal{M}_{1})p(\theta|\phi=\phi_{0}\mathcal{M}_{1})d\theta\\\
&=p(d|\phi=\phi_{0}\mathcal{M}_{1})\,,\end{split}$ (8.36)

in which we used the fact that $\mathcal{M}_{0}$ is nested in
$\mathcal{M}_{1}$ for $\phi=\phi_{0}$. The next step is simply to employ
Bayes’ theorem

$p(\phi=\phi_{0}|d,\mathcal{M}_{1})=\frac{p(d|\phi=\phi_{0}\mathcal{M}_{1})p(\phi=\phi_{0}|\mathcal{M}_{1})}{p(d|\mathcal{M}_{1})}\,,$
(8.37)

and solve for the ratio of model evidences by first substituting equation 8.36
into it to compete the proof.

Practical use of equation 8.32 requires the evaluation the marginal posterior
of $\mathcal{M}_{1}$ at $\phi=\phi_{0}$. As the closed–form expression for the
marginal posterior is rarely available, one can use samples from posterior
(e.g. generated using MCMC) to create a density histogram for $\phi$. Even
better, Kernel Density Estimation (KDE) [158] can be used to approximate the
marginal posterior from samples as,

$\hat{p}_{\mathrm{KDE}}(\phi|d,\mathcal{M}_{1})=\frac{1}{n}\sum_{i=1}^{n}K_{h}\left(\phi-\phi_{i}\right)\,,$
(8.38)

where $K_{h}$ is the kernel and $h$ is the bandwidth, a parameter that
controls the smoothing. The kernel is generally a non–negative function, and
most commonly it is chosen to be a simple Gaussian,

$K_{h}(\phi)=\frac{1}{(2\pi)^{D/2}h^{D}}e^{-\frac{\phi^{2}}{2h^{2}}}\,,$
(8.39)

where $D$ is the dimensionality of $\phi$ (i.e. the number of elements of the
$\phi$ vector). Finally, the value of $h$ can either be determined on the
basis of trial–and–error, or heuristics such as,

$h=1.06\hat{\sigma}n^{-1/5}\,,$ (8.40)

for the $1$–D case where $\hat{\sigma}$ is the standard deviation of the
samples [158].

### 9 Advanced methods

Look on my works, ye Mighty, and despair!
— Percy Shelley, Ozymandias

This chapter introduces two advanced Monte Carlo methods, Sequential Monte
Carlo and Nested sampling, which combine different previously introduced
methods, such as MCMC and importance sampling, in order to provide samples
from posterior distributions and estimate the model evidence. What
distinguishes those two methods from all the previous ones introduced in this
thesis, is their level of complexity and their reliance on multiple individual
algorithms as their constituent parts.

#### 41 Sequential Monte Carlo

Sequential Monte Carlo (SMC) is, from a physics point of view, conceptually
related to the notion of thermodynamic reversibility. For a physical process
starting from a state $A$ and ending in a state $B$, to be thermodynamically
reversible, the transition has to be slow enough such that each intermediate
state of the system is approximately in equilibrium.

##### 41.1 Background

The basic idea of SMC is to slowly guide a population of $n$ particles
$\\{\theta_{i}^{(t)}\\}_{i=1}^{n}$, drawn from a known probability
distribution $\rho(\theta)$, through a series of intermediate distributions
which create a path from $\rho(\theta)$ to the target distribution of interest
$p(\theta)$ [115]. In the context of SMC, the rate of this transition is
governed by the number of intermediate distributions bridging $\rho(\theta)$
to $p(\theta)$. Just like in annealed importance sampling (AIS), SMC relies on
a number of MCMC steps performed in each intermediate step by every particle.
This aims to equilibrate the particles by letting them reach the equilibrium
distribution of each step. Furthermore, when transitioning from an
intermediate distribution to the next, the particle distribution is adjusted
using importance sampling. This guarantees that the particle distribution at
any stage is the correct equilibrium distribution.

The main difference between SMC and AIS is the use of resampling in the case
of SMC. During the run, the particle distribution might experience weight
degeneracy, that is, only a few of the particles have non–negligible
importance weights with the rest of them being vanishingly small. This high
weight–variance can substantially affect any expectation values. In order to
address this issue, SMC performs regular resampling steps, in which the
particle distribution is resampled according to their weights, and the
importance weights are re-initialised to be equal.

SMC methods are particularly suited for challenging target distributions which
exhibit multiple modes. Furthermore, modifications of the main algorithm that
we will present here can also be used for tasks of online learning in which
the data arrive sequentially. These algorithms are most often called by the
name of particle filters [127].

##### 41.2 Bridging the prior and the posterior

A common way to construct such a sequence of intermediate distributions that
bridge a known density $\rho(\theta)$ to the target density $p(\theta)$ is to
interpolate between the two densities

$p_{t}(\theta)\propto\rho^{1-\beta_{t}}(\theta)p^{\beta_{t}}(\theta)\,,\quad
t=1,\dots,m\,,$ (9.1)

where $\beta_{t}$ is a temperature annealing ladder, such that

$0=\beta_{1}<\beta_{2}<\dots<\beta_{m}=1\,.$ (9.2)

In the Bayesian context, a natural choice is to set the prior as the auxiliary
density $\rho(\theta)=p(\theta|\mathcal{M})$ and the posterior as the target
density $p(\theta)=p(\theta|d,\mathcal{M})$. Equation 9.1 then reduces to the
usual annealed or tempered interpolation

$p_{t}(\theta)\propto
p^{\beta_{t}}(d|\theta\mathcal{M})p(\theta|\mathcal{M})\,.$ (9.3)

Although we will focus on this case, the algorithm is valid for any pair of
distributions as long as the support of the auxiliary density encompasses that
of the target.

##### 41.3 Correction – Selection – Mutation

Figure 9.1: Illustration of the Sequential Monte Carlo algorithm with its
three fundamental steps. During the correction step the particles are
reweighted to represent the next probability distribution. Selection removes
the particles with the smaller important weights and multiplies those with
larger weights. Finally, mutation diversifies the particles by moving them.

Given the initial positions of the particles
$\\{\theta_{i}^{(1)}\\}_{i=1}^{n}$ drawn from the prior distribution, as well
as the initial weights $\\{W_{i}^{(1)}\\}_{i=1}^{n}=1/n$, SMC proceeds by the
sequential application of the following three steps, selection, mutation, and
correction, until the posterior density is reached. The procedure that takes
place in a single iteration $t$ is illustrated in Figure 9.1 and involves the
steps:

  1. 1.

Correction / reweighting – During this stage, the weights of the particles are
updated according to

$w_{i}^{(t)}=W_{i}^{(t-1)}\times\frac{p_{t}(\theta_{i}^{(t-1)})}{p_{t-1}(\theta_{i}^{(t-1)})}=W_{i}^{(t-1)}\times\left[p(d|\theta_{i}^{(t-1)},\mathcal{M})\right]^{\beta_{t}-\beta_{t-1}}\,,$
(9.4)

where with $w_{i}$ we denote the unnormalised weights and with
$W_{i}=w_{i}/\sum_{i=1}^{n}w_{i}$ the normalised ones.

The reweighing step accounts and corrects for any deviations of the particle
distribution from the typical set of the target $p_{t}$. The ratio of the
normalisation constants is estimated as

$\frac{\mathcal{Z}_{t}}{\mathcal{Z}_{t-1}}=\sum_{i=1}^{n}w_{i}^{(t)}\,.$ (9.5)

Assuming the density for $t=1$ corresponds to the prior for which
$\mathcal{Z}_{1}=1$, equation 9.5 will eventually lead to the estimation of
the model evidence $\mathcal{Z}_{m}=p(d|\mathcal{M})$.

  2. 2.

Selection / resampling – The particle positions
$\\{\theta_{i}^{(t-1)}\\}_{i=1}^{n}$ are resampled according to their weights
$\\{W_{i}^{(t-1)}\\}_{i=1}^{n}$. The weights are then set again to be equal,
$W_{i}^{(t-1)}=1/n$. Their new, resampled, positions are denoted as
$\\{\tilde{\theta}_{i}^{(t-1)}\\}_{i=1}^{n}$. Particles with small weight
values are removed and those with large importance weights are multiplied.

Resampling can be done using simple multinomial resampling, in which we draw
$n$ new particles, with replacement, with probabilities given by their
weights, or using more advanced schemes characterised by lower variance [113].
This process can be performed in each iteration, or only when some criterion
is triggered (e.g. when the effective sample size of the weights drops below a
threshold).

Finally, caution must be taken when applying resampling too frequently. This
could lead to the phenomenon of weight impoverishment in which there is no
diversity between the particle positions. Fortunately, weight impoverishment
is also reduced by the next step.

  3. 3.

Mutation / propagation – Finally, the population of particles
$\\{\tilde{\theta}_{i}^{(t-1)}\\}_{i=1}^{n}$ is updated and the particles move
to their new positions $\\{\theta_{i}^{(t)}\\}_{i=1}^{n}$ by performing a
number of MCMC steps targeting the density $p_{t}(\theta)$.

The purpose of this step is to diversify the particles and allow their
distribution to approach the stationary distribution. An advantage of SMC is
that the particle distribution from the previous iteration can be used to
construct efficient proposal distributions for MCMC for the current density.
Furthermore, as $n$ particles are updated at once, this step can be done in
parallel. Any MCMC method can be used in this step and there is no requirement
for the final/new positions to be uncorrelated from the initial ones, although
in practice this helps reduce the variance of the estimates.

A common approach is to use the particle covariance $\Sigma^{(t-1)}$ to
construct a normal proposal distribution
$q(\theta^{\prime}|\theta)=\mathcal{N}(\theta^{\prime}|\theta,\Sigma^{(t-1)})$.

Once all three steps are completed, the value of $\beta$ is updated and the
process is repeated again until $\beta$ reaches the value of one. The names of
those three steps are inspired by natural selection and evolutionary
programming. The reason is the apparent analogy with genetic algorithms [106].
More specifically, reweighting, resampling and propagation have the roles of
correction, selection and mutation in genetic algorithms, in which the
particle positions are the genes and the importance weights play the role of
the so–called fitness. A critical difference with most genetic algorithms is
the fact that SMC solves a sampling task, not an optimisation one, and thus
the solution is represented by the distribution of the particles and not by
any particle individually.

Algorithm 13 Sequential Monte Carlo

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})$, and a local MCMC kernel
$\theta^{\prime}\leftarrow\mathcal{K}(\theta;f(\theta))$ (e.g. $N$ steps of
random–walk Metropolis update)

0: Posterior samples and estimate of the model evidence $\mathcal{Z}$

1: Initialise temperature parameter $\beta_{1}=1$

2: Initialise estimate of evidence $\mathcal{Z}=1$

3: for $i=1$ to $n$ do

4: Draw particle positions from the prior $\theta_{i}^{(1)}\sim\pi(\theta)$

5: Initialise particle weights $W_{i}^{(1)}=1/n$

6: end for

7: while $\beta_{t}\neq 1$ do

8: Update iteration index $t\leftarrow t+1$

9: Set temperature $\beta_{t}$ solving
$\left(\sum_{i=1}^{n}w_{i}^{(t)}(\beta_{t})\right)^{2}/\sum_{i=1}^{n}\left(w_{i}^{(t)}(\beta_{t})\right)^{2}=\alpha\times
n$ where the importance weights are computed as $w_{i}^{(t)}\leftarrow
W_{i}^{(t-1)}\mathcal{L}(\theta_{i}^{(t-1)})^{\beta_{t}-\beta_{t-1}}$

10: Update evidence estimate $\mathcal{Z}\leftarrow\mathcal{Z}\times
n^{-1}\sum_{i=1}^{n}w_{i}^{(t)}$

11: $\left\\{\tilde{\theta}_{i}^{(t-1)}\right\\}_{i=1}^{n}\leftarrow$ resample
$\left\\{\theta_{i}^{(t-1)}\right\\}_{i=1}^{n}$ according to
$\left\\{w_{i}^{(t)}\right\\}_{i=1}^{n}$

12: for $i=1$ to $n$ do

13: Reset weights $W_{i}^{(t)}\leftarrow 1/n$

14: end for

15: Update particles using MCMC
$\left\\{\theta_{i}^{(t)}\right\\}_{i=1}^{n}\leftarrow\mathcal{K}\left(\left\\{\tilde{\theta}_{i}^{(t-1)}\right\\}_{i=1}^{n};\pi(\theta)\mathcal{L}(\theta)^{\beta_{t}}\right)$

16: end while

##### 41.4 Effective sample size

A common measure of the quality of the importance weights of the particles, at
any iteration of the SMC run, is the effective sample size (ESS)

$\text{ESS}_{t}=\frac{\mathbb{E}_{p_{t}}\left[w^{(t)}\right]^{2}}{\mathbb{E}_{p_{t}}\left[\left(w^{(t)}\right)^{2}\right]}\,.$
(9.6)

which can be estimated as:

$\hat{\text{ESS}}_{t}=\frac{\left(\sum_{i=1}^{n}w_{i}^{(t)}\right)^{2}}{\sum_{i=1}^{n}(w_{i}^{(t)})^{2}}=\frac{1}{\sum_{i=1}^{n}(W_{i}^{(t)})^{2}}\,.$
(9.7)

##### 41.5 Setting the temperature ladder

A $\beta$ ladder can be specified a priori or determined adaptively during the
run [61]. In the first case, the resampling step is usually triggered whenever
the ESS drops below a prespecified threshold value (e.g. $50\%-95\%$). In the
latter case, the next value of $\beta$ is chosen adaptively such that the ESS
has an approximately constant fraction $\alpha$ (e.g. $50\%-95\%$) of the
number of particles $n$ throughout the duration of the SMC run. Numerically,
this can be done by solving

$\frac{\left[\sum_{i=1}^{n}w_{i}^{(t)}(\beta_{t})\right]^{2}}{\sum_{i=1}^{n}\left[w_{i}^{(t)}(\beta_{t})\right]^{2}}=\alpha\times
n\,,$ (9.8)

for the next $\beta_{t}$ such that $\beta_{t-1}<\beta_{t}\leq 1$ using, for
instance, the _bisection method_ [23].

#### 42 Nested sampling

Nested sampling (NS), originally developed by [160, 161], is a method for
estimating the model evidence $\mathcal{Z}=p(d|\mathcal{M})$. The basic idea
is to approximate the evidence by integrating the prior in nested shells of
constant likelihood. Despite its original purpose to estimate the model
evidence, NS can also provide weighted samples from the posterior distribution
as an optional byproduct. Therefore, the method is suitable for both tasks of
parameter estimation and model comparison [4]. Over the years, many variants
of NS have emerged, with each one aiming to improve a different aspect of the
original version [20, 44, 81].

##### 42.1 Multi–dimensional integration

Figure 9.2: Illustration comparing two ways which one can use to approximate
the model evidence integral. The left panel shows the direct multi–dimensional
integration over the parameters. The right panel shows the one–dimensional
integration over the prior volume $X$ enclosed in the iso–likelihood contours.

NS attempts to compute the evidence integral,

$\mathcal{Z}=\int\mathcal{L}(\theta)\pi(\theta)d\theta\,,$ (9.9)

where $\mathcal{L}(\theta)=p(d|\theta,\mathcal{M})$ is the likelihood function
and $\pi(\theta)=p(\theta|\mathcal{M})$ is the prior, by transforming it into
a one–dimensional integral over the prior volume

$X(\lambda)=\int_{\mathcal{L}(\theta)>\lambda}\pi(\theta)d\theta,$ (9.10)

enclosed in the iso–likelihood contour defined by
$\mathcal{L}(\theta)=\lambda$. Equation 9.9 can then be written as

$\mathcal{Z}=\int_{0}^{+\infty}X(\lambda)d\lambda=\int_{0}^{1}\mathcal{L}(X)dX\,,$
(9.11)

assuming that $\mathcal{L}(X(\lambda))=\lambda$ exists. Figure 9.2 illustrates
the equivalency between the integrals of equations 9.9 and 9.11. Unlike
equation 9.9, the above integral is now 1–dimensional and can be approximated
using standard numerical integration techniques (e.g. quadrature),

$\mathcal{Z}=\sum_{i=1}^{m}\mathcal{L}_{i}w_{i}\,,$ (9.12)

where,

$w_{i}=\frac{X_{i-1}-X_{i+1}}{2}\,.$ (9.13)

This of course assumes that we are able to evaluate the iso–likelihood
contours $\mathcal{L}_{i}=\mathcal{L}(X_{i})$ associated with an ordered
collection of samples with prior volume $1>X_{1}>X_{2}>\dots>X_{m}>0$. This is
illustrated in Figure 9.3 for a collection of $8$ samples that are uniformly
distributed in the prior volume.

Figure 9.3: Illustration of $8$ samples drawn uniformly from the prior with
their respective iso–likelihood contours (left), along with their
corresponding contributions to the evidence integral (right).

Using the simpler weights $w_{i}=X_{i}-X_{i+1}$ in equation 9.12 a lower bound
on the evidence can be estimated as

$\mathcal{Z}\geq\sum_{i=1}^{m}\mathcal{L}_{i}\left(X_{i}-X_{i+1}\right)\,.$
(9.14)

Similarly, an upper bound also exists, using $w_{i}=X_{i-1}-X_{i}$, which can
be written as

$\mathcal{Z}\leq\sum_{i=1}^{m}\mathcal{L}_{i}\,\left(X_{i-1}-X_{i}\right)+X_{m}\mathcal{L}_{\mathrm{max}}\,,$
(9.15)

where $\mathcal{L}_{\mathrm{max}}$ is the maximum likelihood value to be found
as $X\to 0$.

Soon after its original inception, it was realised that a NS run can also be
used for the task of parameter estimation without any additional computation.
In particular, the collected samples combined with their normalised weights

$p_{i}=\frac{\mathcal{L}_{i}w_{i}}{\mathcal{Z}}\,,$ (9.16)

correspond to weighted samples from the posterior distribution and thus can be
used to compute expectation values

$\mathbb{E}_{p}\left[f(\theta)\right]=\sum_{i}^{n}p_{i}f(\theta_{i})\,.$
(9.17)

##### 42.2 Sampling procedure

Figure 9.4: Illustration of the nested sampling procedure. Given some
uniformly distributed points from the prior, we identify and remove the worst
point, that is, the point with the minimum likelihood value.
$\mathcal{L}_{min}$, and replace it a new point sampled from the prior subject
to the likelihood constrain $\mathcal{L}>\mathcal{L}_{min}$. Finally, the
volume is contracted to account for the removal of the worst point.

The NS algorithm begins by drawing a collection of points
$\\{\theta_{i}\\}_{i=1}^{n}$ uniformly from the prior, often called live
points. We can associate each live point $\theta_{i}$ with a prior volume
$X$–value, namely the volume that would be enclosed by the iso–likelihood
contour $\mathcal{L}_{i}=\mathcal{L}(\theta_{i})$. On average, we expect
roughly half of the live points to fall inside the iso–likelihood contour
corresponding to half prior volume $X=1/2$, one quarter to $X=1/4$, one eighth
to $X=1/8$ and so on. In other words, since the live points are uniformly
distributed under the prior, the corresponding $X$–values are uniformly
distributed between $0$ and $1$. This is illustrated in Figure 9.3 and the top
panel of Figure 9.4.

What we described so far is only the first step of the algorithm, and one
still needs a way to propagate the live points into regions of smaller prior
volume (i.e. lower $X$) in order to probe iso–likelihood contours
corresponding to higher likelihood values. NS achieves this by first
identifying the live point with the lowest likelihood value
$\mathcal{L}^{*}=\mathcal{L}_{1}$, corresponding to volume $X_{1}$ and
removing it. The remaining live points are now distributed over a compressed
volume $X_{1}$. On average, the volume compression factor is

$t=e^{-1/n}\,,$ (9.18)

such that the compressed volume is $X_{1}=t\times X_{0}$, where $X_{0}=1$ is
the initial total volume. Finally, we sample a new live point to replace the
one that we removed. The new point is sampled uniformly from the prior subject
to the constrain $\mathcal{L}>\mathcal{L}^{*}$, that is, from the
likelihood–constrained prior

$\pi^{*}(\theta)\propto\begin{cases}\pi(\theta)&\mathrm{if}\mathcal{L}(\theta)>\mathcal{L}^{*}\,,\\\
0&\text{otherwise}\,.\end{cases}$ (9.19)

This whole process, that is shown in Figure 9.4, is repeated multiple times
until a criterion for termination is met. In each iteration, the volume
shrinks on average by the compression factor of equation 9.18.

##### 42.3 Termination criterion

During an NS run, the remaining prior volume $X$ asymptotically approaches
$0$. The fact that we can only perform a finite number of steps means that it
is unavoidable to introduce a truncation error into the evidence estimate of
equation 9.11. A common way of determining when to stop is to approximately
estimate the amount of remaining evidence and terminate the run when this can
be considered negligible for the purpose of the analysis.

Perhaps the simplest way to roughly estimate the remaining evidence is by
utilising the upper bound of equation 9.15. In this case, the remaining
evidence is approximated as
$\Delta\mathcal{Z}=\mathcal{L}_{\mathrm{max}}X_{m}$, where
$\mathcal{L}_{\mathrm{max}}$ is the maximum likelihood value, estimated from
the current population of live points, and $X_{m}$ is simply the estimate of
the remaining prior volume. An alternative would be to use the mean likelihood
of the live points and get $\Delta\mathcal{Z}=\bar{\mathcal{L}}X_{m}$. The run
then terminates when $\Delta\mathcal{Z}/\mathcal{Z}$ drops below a
prespecified threshold.

Of course, neither of these approaches guarantees that the run will terminate
early and that beyond lies a “spike” of huge likelihood. Upon deciding to
stop, however, the current estimate of the model evidence is approximately
corrected by either adding $\Delta\mathcal{Z}$ or removing the live points
one–by–one in accordance with the NS procedure and adding their respective
evidence contributions $\mathcal{L}_{i}X_{i}$, but without replacing them with
new ones.

##### 42.4 Uncertainty

So far we have assumed that the compression factor is given by equation 9.18,
however, this is simply the mean compression factor associated with the
removal of the outermost or lowest–likelihood live point. In truth, the prior
volume bounded by the iso–likelihood contour of that point can be slightly
different from what the mean compression factor predicts. The compression in
volume $t$ associated with the removal of the outermost of $n$ live points
follows a $\mathrm{Beta}(n,1)$ probability distribution with density,

$p(t)=nt^{n-1}\,,$ (9.20)

where the first $n$ factor comes from the fact that any live point could be
the outermost, and the second factor from the fact that the remaining $n-1$
live points lie uniformly distributed above the outermost at $t$.

The compression factors can therefore be sampled from the probability
distribution of equation 9.20 instead of just assuming their expected value of
equation 9.18. Furthermore, we can use equation 9.20 to compute the
expectation values

$\mathbb{E}[\log t]=-\frac{1}{n}\,,\quad\,\mathrm{Var}[\log
t]=\frac{1}{n^{2}}\,.$ (9.21)

Since the individual $\log t$ are independent, we expect that after $i$ steps,
the prior volume to have shrunk to

$\log X_{i}=-\frac{i\pm\sqrt{i}}{n}\,.$ (9.22)

What the above expression means is that there is uncertainty in the estimates
of the compression factor which enter into the prior volume estimates too. In
other words, there is uncertainty in the number of steps $i$ required for the
prior volume to shrink to a certain value $X_{i}$.

The uncertainty originating from the noisy estimates of the compression factor
$t$ also propagates into the estimate of the model evidence. To quantify this
we need to consider the information gained when transitioning from the prior
$\pi(\theta)$ to the posterior $p(\theta)$, given by the Kullback–Leibler (KL)
divergence,

$H=\int p(\theta)\log\left(\frac{p(\theta)}{\pi(\theta)}\right)d\theta\,.$
(9.23)

We can rewrite the above equation in terms of the prior volume $X$, as

$\begin{split}H&=\int p(X)\log p(X)dX\\\ &=-\int p(X)\log XdX+\int p(\log
X)\log p(\log X)d\log X\,,\end{split}$ (9.24)

where $p(X)=\mathcal{L}(X)/\mathcal{Z}$ is the volume posterior density.
Ignoring the second term on the right hand size, which is subdominant, we thus
get that the KL divergence provides a measure of the compression we require to
reach the bulk of the posterior mass,

$H=-\log X\,.$ (9.25)

Comparing equation 9.22 and 9.25 we roughly expect $nH\pm\sqrt{nH}$ steps to
reach the bulk of the posterior mass. Equivalently, the uncertainty introduced
into the estimate of the model evidence is

$\Delta\log\mathcal{Z}=\sqrt{\frac{H}{n}}\,.$ (9.26)

Of course, the above expression does not include any sources of numerical
error such as truncation error.

##### 42.5 Likelihood–constrained prior sampling

The efficient application of the NS algorithm requires sampling from the prior
distribution $\pi(\theta)$ subject to the likelihood constrain
$\mathcal{L}>\mathcal{L}^{*}$. Unfortunately, drawing points from the prior
until the likelihood criterion is met is not feasible in practice, as the
volume contained in the constrained prior shrinks exponentially with each
iteration. For this reason, two different approaches, region and step samplers
are often employed in order to produce samples from the likelihood–constrained
prior.

For the sake of simplicity, both samplers usually operate in the latent
parameter space that the prior is uniform over the unit hypercube. In this
case, the practitioner specifies their prior preference, not by providing a
(log–) probability density function, but by defining the inverse–cumulative
density function $\Phi^{-1}$ that transforms points $\phi$ in the hypercube to
points $\theta$ in the original parameter space. For instance, let us assume
that we require a normal prior on a parameter
$\theta\sim\mathcal{N}(\mu,\sigma^{2})$. We can transform a unit hypercube
parameter $\phi\sim\mathcal{U}(0,1)$, using the standard normal distribution’s
inverse–cumulative density function $\Phi^{-1}$, such that,

$\theta=\mu+\Phi^{-1}(\phi)\sigma\,.$ (9.27)

###### 42.5.1 Region samplers

The basic idea behind region samplers is to construct a hypersurface that
bounds a given iso–likelihood contour. In practice, this is done using simple
geometric shapes (e.g. spheres, ellipses, etc.). The hypersurface must
encompass the current distribution of live points and at least contain the
currently estimated volume. One can then sample uniformly from within the
volume enclosed by the hypersurface until the likelihood–constraint is
satisfied. In order to reduce the risk of missing parts of the currently
estimated volume, the bounding region is usually expanded by a prespecified
factor or using cross–validation of the live points.

Most region samplers attempt to construct such a bounding region by wrapping
the current generation of live points with one or multiple ellipsoids. Using
multiple ellipsoids offers some flexibility in the case of multimodal
posterior distributions. The most popular such sampler is the MultiNest
algorithm that determines the shape and location of the ellipsoids based on
the mean and covariance of the live points, by first estimating the number of
distinct modes, and thus required ellipsoids, using a clustering algorithm
[46].

Region samplers have to face serious challenges when the complexity of the
posterior or the dimensionality of the parameter space increases. In the first
case, the ability to accurately bound the current volume depends on the number
of live points, with a higher number often resulting in better bounding
regions. The second limitation arises from the curse of dimensionality. As the
number of dimensions increases, most of the volume of the bounding shape
concentrates near its edges, and given that the bounding region is often
chosen to be significantly broader than the encompassing
likelihood–constrained volume to guarantee that no parts are encroached, the
total number of samples until one is found to lie within the iso–likelihood
contour increases exponentially. As a consequence, region samplers are more
efficient and appropriate for low–dimensional problems, that is, $D\lesssim
10-20$.

###### 42.5.2 Step samplers

On the other hand, step samplers do not rely on a bounding region and thus
bypass some of the pathologies of region samplers. Instead, they evolve a
randomly chosen live point through a sequence of local steps to an
approximately independent position. This is usually achieved using some MCMC
method targeting the likelihood–constrained prior of equation 9.19 as the
target distribution. The advantage of using MCMC methods in the context of NS
is that, in each iteration, one can use the distribution of the live points to
construct effective proposal distributions for the MCMC sampler.

Although step samplers enjoy a better scaling with the number of dimensions
than region samplers, there are still challenges in their use. First of all,
determining the minimum number of steps to perform for the new point to be
independent of its starting position (i.e. the randomly chosen live point) is
not trivial. Although small correlations can be effectively ignored, larger
violations can have catastrophic results and lead to substantial bias in the
final estimates of NS. Furthermore, the step sampler must be tuned
appropriately to achieve good sampling performance. Adaptation during a given
iteration has to be diminishing in order to avoid spurious effects and biases.

Although any MCMC method (e.g. random walk Metropolis, slice sampling, etc.)
can in principle be used as a step sampler, there are also methods that are
naturally suited and have been developed for use in the context of sampling
from the likelihood–constrained prior. One such example is Galilean Monte
Carlo (GMC) [47, 162, 163] that samples by moving consistently along a
direction until a proposed point is rejected, by being outside the
iso–likelihood contour. In this case, the sampler reflects off the current
iso–likelihood boundary.

##### 42.6 Parallelisation

Parallelising NS is not as straightforward as with other Monte Carlo
algorithms (e.g. Sequential Monte Carlo), as the method relies on updating a
single (worst) point at a time. In general, we would like to generate as many
candidate points per step as the number of available CPUs (i.e.
$n_{\mathrm{CPU}}$) and evaluate their likelihoods in parallel. In this case,
there are three strategies that one can follow:

  1. 1.

Replace a single live point and discard as many as $n_{\mathrm{CPU}}-1$
acceptable live points. This scheme is quite wasteful, particularly in cases
in which it is likely that more than one candidate point satisfies the
likelihood constraint.

  2. 2.

Replace the worst (i.e. lowest likelihood) $n_{\mathrm{CPU}}$ live points in a
single step. This results in linear speed–up with respect to the number of
CPUs but increases the variance of the evidence estimate by a factor of
$\sqrt{n_{\mathrm{CPU}}}$.

  3. 3.

Replace a single live point and consider the other $n_{\mathrm{CPU}}-1$
candidates for subsequent iterations. This results in a speed–up of
$n\log(1+n_{\mathrm{CPU}}/n)$ which is approximately linear for
$n_{\mathrm{CPU}}<<n$. The “diminishing returns” represented by the
logarithmic factor in this expression originate from the fact that the
likelihood threshold increases as the run progresses, and thus the points
might not be valid for a subsequent iteration. This strategy is the most
widely employed in practice.

Finally, it is important to note that apart from parallelising a single NS
run, it is also possible to combine different, possibly parallel, independent
NS runs into a joint one, thus achieving linear scaling. In order to combine
two or more runs together, we collect the points from all runs as live points
and begin by removing the worst point, which with no loss of generality we
assume that it belongs to run A. Then, as a replacement that satisfies the
likelihood constraint, we simply take the replacement that was originally used
in run A. We then proceed with the next worst point until all points are
accounted for.

Algorithm 14 Nested sampling

0: termination criterion (e.g. $\Delta\log\mathcal{Z}\leq\epsilon)$, number of
live points $n$, an estimate of the compression factor e.g. $t=\exp(-1/n)$,
prior distribution $\pi(\theta)$, likelihood function
$p_{i}=\mathcal{L}(\theta)$)

0: Estimate of model evidence $\mathcal{Z}$, posterior samples
$\tilde{\theta}_{i}$ with weights $\mathcal{L}_{i}w_{i}/\mathcal{Z}$

1: Initialise volume $X=1$ and evidence $\mathcal{Z}=0$

2: Draw $n$ live points from the prior
$\theta_{1},\theta_{2},\dots,\theta_{n}\sim\pi(\theta)$

3: repeat

4: Find the minimum likelihood value of the live points
$\mathcal{L}^{*}\leftarrow\min\left(\mathcal{L}(\theta_{1}),\dots,\mathcal{L}(\theta_{n})\right)$

5: Replace live point $\theta^{*}$ corresponding to $\mathcal{L}^{*}$ with a
new point from the prior satisfying $\mathcal{L}>\mathcal{L}^{*}$

6: Set $w^{*}\leftarrow\Delta X$ where $\Delta X=(1-t)X$

7: Update estimate of the evidence
$\mathcal{Z}\leftarrow\mathcal{Z}+\mathcal{L}^{*}w^{*}$

8: Store values of $w_{i}\leftarrow w^{*}$,
$\mathcal{L}_{i}\leftarrow\mathcal{L}^{*}$, and
$\tilde{\theta}_{i}\leftarrow\theta^{*}$

9: Contract volume $X\leftarrow tX$

10: until termination criteria satisfied

11: Update evidence
$\mathcal{Z}\leftarrow\mathcal{Z}+\frac{1}{n}\sum_{j=1}^{n}\mathcal{L}(\theta_{j})\,X$

##  Part III Novel Developments

### 10 Ensemble Slice Sampling

This chapter presents Ensemble Slice Sampling which is the main contribution
introduced in the paper titled Ensemble Slice Sampling: Parallel, black–box,
and gradient–free inference that was published in the journal Statistics and
Computing in 2021 [95]. The content of the chapter is almost identical to that
included in the aforementioned publication with the exception of minor text
and figure formatting differences.

Slice Sampling has emerged as a powerful Markov Chain Monte Carlo algorithm
that adapts to the characteristics of the target distribution with minimal
hand-tuning. However, Slice Sampling’s performance is highly sensitive to the
user-specified initial length scale hyperparameter and the method generally
struggles with poorly scaled or strongly correlated distributions. This paper
introduces Ensemble Slice Sampling (ESS), a new class of algorithms that
bypasses such difficulties by adaptively tuning the initial length scale and
utilising an ensemble of parallel walkers in order to efficiently handle
strong correlations between parameters. These affine–invariant algorithms are
trivial to construct, require no hand-tuning, and can easily be implemented in
parallel computing environments. Empirical tests show that Ensemble Slice
Sampling can improve efficiency by more than an order of magnitude compared to
conventional MCMC methods on a broad range of highly correlated target
distributions. In cases of strongly multimodal target distributions, Ensemble
Slice Sampling can sample efficiently even in high dimensions. We argue that
the parallel, black-box and gradient-free nature of the method renders it
ideal for use in scientific fields such as physics, astrophysics and cosmology
which are dominated by a wide variety of computationally expensive and non-
differentiable models.

#### 43 Introduction

_Bayesian inference and data analysis_ has become an integral part of modern
science. This is partly due to the ability of Markov Chain Monte Carlo (MCMC)
algorithms to generate samples from intractable probability distributions.
MCMC methods produce a sequence of samples, called a _Markov chain_ , that has
the target distribution as its equilibrium distribution. The more samples are
included, the more closely the distribution of the samples approaches the
target distribution. The Markov chain can then be used to numerically
approximate expectation values (e.g. parameter uncertainties, marginalised
distributions).

Common MCMC methods entail a significant amount of time spent hand-tuning the
hyperparameters of the algorithm to optimize its efficiency with respect to a
target distribution. The emerging and routine use of such mathematical tools
in science calls for the development of black-box MCMC algorithms that require
no hand-tuning at all. This need led to the development of adaptive MCMC
methods like the Adaptive Metropolis algorithm [77] which tunes its proposal
scale based on the sample covariance matrix. Unfortunately, most of those
algorithms still include a significant number of hyperparameters (e.g.
components of the covariance matrix) rendering the adaptation noisy.
Furthermore, the tuning is usually performed on the basis of prior knowledge,
such as one or more long preliminary runs which further slow down the
sampling. Last but not least, there is no reason to believe that a single
Metropolis proposal scale is optimal for the whole distribution (i.e. the
appropriate scale could vary from one part of the distribution to another).
Another approach to deal with those issues would be to develop methods that by
construction require no or minimal hand-tuning. An archetypal such method is
the Slice Sampler [133], which has only one hyperparameter, the initial length
scale.

It should be noted that powerful adaptive methods that require no hand-tuning
(although they do require preliminary runs) already exist. Most notable of
them is the No U-Turn Sampler (NUTS) [84], an adaptive extension of
Hamiltonian Monte Carlo (HMC) [135]. However, such methods rely on the
gradient of the log probability density function. This requirement is the
reason why these methods are limited in their application in quantitative
fields such as physics, astrophysics and cosmology, which are dominated by
computationally costly non–differentiable models. Thus, our objective in this
paper is to introduce a parallel, black-box and gradient–free method that can
be used in the aforementioned scientific fields.

This paper presents Ensemble Slice Sampling (ESS), an extension of the
Standard Slice Sampling method. ESS naturally inherits most of the benefits of
Standard Slice Sampling, such as the acceptance rate of $1$, and most
importantly the ability to adapt to the characteristics of a target
distribution without any hand-tuning at all. Furthermore, we will show that
ESS’s performance is insensitive to linear correlations between the
parameters, thus enabling efficient sampling even in highly demanding
scenarios. We will also demonstrate ESS’s performance in strongly multimodal
target distributions and show that the method samples efficiently even in high
dimensions. Finally, the method can easily be implemented in parallel taking
advantage of multiple CPUs thus facilitating Bayesian inference in cases of
computationally expensive models.

Our implementation of ESS is inspired by [174]. However, our method improves
upon that by extending the direction choices (e.g. Gaussian and global move),
adaptively tuning the initial proposal scale, and parallelising the algorithm.
[138] developed a general algorithm based on the elliptical slice sampling
method [126] and a Gaussian mixture approximation to the target distribution.
ESS utilises an ensemble of parallel and interacting chains, called walkers.
Other methods that are based on the ensemble paradigm include the Affine
Invariant Ensemble Sampler [72] and the Differential Evolution MCMC [169]
along with its various extensions [170, 182], as well as more recent
approaches that are based on langevin diffusion dynamics [53, 54] and the time
discretization of stochastic differential equations [110] in order to achieve
substantial speedups.

In Section 44, we will briefly discuss the Standard Slice Sampling algorithm.
In Section 45, we will introduce the Ensemble Slice Sampling method. In
Section 46 we will investigate the empirical evaluation of the algorithm. We
reserve Sections 47 and 48 for discussion and conclusion, respectively.

#### 44 Standard Slice Sampling

Slice Sampling is based on the idea that sampling from a distribution $p(x)$
whose density is proportional to $f(x)$ is equivalent to uniformly sampling
from the region underneath the graph of $f(x)$. More formally, in the
univariate case, we introduce an auxiliary variable, the height $y$, thus
defining the joint distribution $p(x,y)$, which is uniform over the region
$U=\\{(x,y):0<y<f(x)\\}$. To sample from the marginal density for $x$, $p(x)$,
we sample from $p(x,y)$ and then we ignore the $y$ values.

Generating samples from $p(x,y)$ is not trivial, so we might consider defining
a Markov chain that will converge to that distribution. The simplest, in
principle, way to construct such a Markov chain is via Gibbs sampling. Given
the current $x$, we sample $y$ from the conditional distribution of $y$ given
$x$, which is uniform over the range $(0,f(x))$. Then we sample the new $x$
from the slice $S=\\{x:y<f(x)\\}$.

Generating a sample from the slice $S$ may still be difficult, since we
generally do not know the exact form of $S$. In that case, we can update $x$
based on a procedure that leaves the uniform distribution of $S$ invariant.
[133] proposed the following method:

Given the current state $x_{0}$, the next one is generated as:

  1. 1.

Draw $y_{0}$ uniformly from $(0,f(x_{0}))$, thus defining the horizontal slice
$S=\\{x:y_{0}<f(x)\\}$,

  2. 2.

Find an interval $I=(L,R)$ that contains all, or much, of $S$ (e.g. using the
stepping-out procedure defined below),

  3. 3.

Draw the new point $x_{1}$ uniformly from $I\cap S$.

Figure 10.1: The plot shows the univariate slice sampling method. Given an
initial value $x_{0}$, a value $y_{0}$ is uniformly sampled along the vertical
slice $(0,f(x_{0}))$ (green dashed line) thus defining the initial point (blue
star). An interval $(L,R)$ is randomly positioned horizontally around the
initial point, and then it is expanded in steps of size $\mu=R-L$ until both
of its ends $L^{\prime},R^{\prime}$ are outside the slice. The new point
(green star) is generated by repeatedly sampling uniformly from the expanded
interval $(L^{\prime},R^{\prime})$ until a point is found inside the slice.
Points outside the slice (e.g. the red star) are used to shrink the interval
$(L^{\prime},R^{\prime})$ by moving $L^{\prime}$ or in this case $R^{\prime}$
to that point and accelerate the sampling procedure.

In order to find the interval $I$, [133] proposed to use the _stepping-out_
procedure that works by randomly positioning an interval of length $\mu$
around the point $x_{0}$ and then expanding it in steps of size $\mu$ until
both ends are outside of the slice. The new point $x_{1}$ is found using the
_shrinking_ procedure, in which points are uniformly sampled from $I$ until a
point inside $S$ is found. Points outside $S$ are used to shrink the interval
$I$. The stepping-out and shrinking procedures are illustrated in Figure 10.1.
By construction, the stepping-out and shrinking procedures can adaptively tune
a poor estimate of the length scale $\mu$ of the initial interval. The length
scale $\mu$ is the only free hyperparameter of the algorithm. For a detailed
review of the method we direct the reader to [133] and [117] (also Exercise
30.12 in that text).

It is important to mention here that for multimodal distributions there is no
guarantee that the slice would cross any of the other modes, especially if the
length scale is underestimated initially. Ideally, in order to provide a large
enough initial value of the scale factor $\mu$, prior knowledge of the
distance between the modes is required. As we will show in the next section,
Ensemble Slice Sampling does not suffer from this complication and can handle
strongly multimodal distributions efficiently.

#### 45 Ensemble Slice Sampling

The univariate slice sampling scheme can be used to sample from multivariate
distributions by sampling repeatedly along each coordinate axis in turn (one
parameter at a time) or by sampling along randomly selected directions [117].
Using either of those choices, the Standard Slice Sampler performs acceptably
in cases with no strong correlations in parameter space. The overall
performance of the algorithm generally depends on the number of expansions and
contractions during the stepping-out and shrinking procedures, respectively.
Ideally we would like to minimize that number. A reasonable initial estimate
of the length scale is still required in order to reduce the amount of time
spent expanding or contracting the initial interval.

However, when strong correlations are present two issues arise. First, there
is no single value of the initial length scale that minimizes the
computational cost of the stepping-out and shrinking procedures along all
directions in parameter space. The second problem concerns the choice of
direction. In particular, neither the component-wise choice (one parameter at
a time) nor the random choice is suitable in strongly correlated cases. Using
such choices results in highly autocorrelated samples.

Our approach would be to target each of those two issues individually. The
resulting algorithm, Ensemble Slice Sampling (ESS), is invariant under affine
transformations of the parameter space, meaning that its performance is not
sensitive to linear correlations. Furthermore, ESS minimizes the computational
cost of finding the slice by adaptively tuning the initial length scale. Last
but not least, unlike most MCMC methods, ESS is trivially parallelizable, thus
enabling the data analyst to take advantage of modern high performance
computing facilities with multiple CPUs.

##### 45.1 Adaptively tuning the length scale

Let us first consider the effect of the initial length scale on the
performance of the univariate slice sampling method. For instance, if the
initial length scale is $\lambda$ times smaller than the actual size of the
slice, then the stepping-out procedure would require $\mathcal{O}(\lambda)$
steps in order to fix this. However, in this case, since the final interval is
an accurate approximation of the slice there would probably be no contractions
during the shrinking phase. On the other hand, when the initial length scale
is larger than the actual slice then the number of expansions would be either
one or zero. In this case though, there would be a number of contractions.

###### 45.1.1 Stochastic approximation

As the task is to minimize the total number of expansions and contractions we
employ and adapt the _Robbins–Monro_ stochastic approximation algorithm [148]
of [171]. Ideally, based on the reasoning of the previous paragraph, only one
expansion and one contraction will take place. Therefore, the target ratio of
number of expansions to total number of expansions and contractions is $1/2$.
To achieve this, we update the length scale $\mu$ based on the following
recursive formula:

$\mu^{(t+1)}=2\mu^{(t)}\frac{N_{e}^{(t)}}{N_{e}^{(t)}+N_{c}^{(t)}}\,,$ (10.1)

where $N_{e}^{(t)}$ and $N_{c}^{(t)}$ are the number of expansions and
contractions during iteration $t$. It is easy to see that when the fraction
$N_{e}^{(t)}/(N_{e}^{(t)}+N_{c}^{(t)})$ is larger than $1/2$ the length scale
$\mu$ will be increased. In the case where the fraction is smaller than $1/2$
the length scale $\mu$ will be decreased accordingly. The optimization can
stop either when the fraction is close to $1/2$ within a threshold or when a
maximum number of tuning steps has been completed. The pseudocode for the
first case is shown in Algorithm 15. In order to preserve detailed balance it
is important to be sure that the adaptation stops after a finite number of
iterations. In practice this happens after $\mathcal{O}(10)$ iterations. An
alternative would be to use diminishing adaptation [152] but we found that our
method is sufficient in practice (see Section 4.3 for more details).

Algorithm 15 Function to tune the length scale $\mu$.

1: function TuneLengthScale($t$, $\mu^{(t)}$, $N_{e}^{(t)}$, $N_{c}^{(t)}$,
$M^{\text{adapt}}$)

2: if $t\leq M^{\text{adapt}}$ then

3: Compute $\mu^{(t+1)}$ using Equation 10.1,

4: return $\mu^{(t+1)}$

5: else

6: return $\mu^{(t)}$

7: end if

##### 45.2 The choice of direction and parallelisation

In cases where the parameters are correlated we can accelerate mixing by
moving more frequently along certain directions in parameter space. One way of
achieving this is to exploit some prior knowledge about the covariance of the
target distribution. However, such an approach would either require
significant hand-tuning or noisy estimations of the sample covariance matrix
during an initial run of the sampler. For that reason we employ a different
approach to exploit the covariance structure of the target distribution and
preserve the hand-tuning-free nature of the algorithm.

###### 45.2.1 Ensemble of walkers

Following the example of [72] we define an ensemble of parallel chains, called
walkers. In our case though, each walker is an individual slice sampler. The
sampling proceeds by moving one walker at a time by slice sampling along a
direction defined by a subset of the rest of the walkers of the ensemble. As
long as the aforementioned direction does not depend on the position of the
current walker, the resulting algorithm preserves the detailed balance of the
chain. Moreover, assuming that the distribution of the walkers resembles the
correlated target distribution, the chosen direction will prefer directions of
correlated parameters.

We define an ensemble of $N$ parallel walkers as the collection
$S=\\{\mathbf{X_{1}},\dots,\mathbf{X_{N}}\\}$. The position of each individual
walker $\mathbf{X_{k}}$ is a vector $\mathbf{X_{k}}\in\mathbb{R}^{D}$ and
therefore we can think of the ensemble $S$ as being in $\mathbb{R}^{ND}$.
Assuming that each walker is drawn independently from the target distribution
with density $p$, then the target distribution for the ensemble would be the
product

$P(\mathbf{X_{1}},\dots,\mathbf{X_{N}})=\prod_{k=1}^{N}p(\mathbf{X_{k}})\,.$
(10.2)

The Markov chain of the ensemble would preserve the product density of
equation 10.2 without the individual walker trajectories being Markov. Indeed,
the position of $\mathbf{X_{k}}$ at iteration $t+1$ can depend on
$\mathbf{X_{j}}$ at iteration $t$ with $j\neq k$.

Given the walker $\mathbf{X_{k}}$ that is to be updated there are arbitrary
many ways to define a direction vector from the complementary ensemble
$S_{[k]}=\\{\mathbf{X_{j}},\>\forall j\neq k\\}$. Here we will discuss a few
of them. Following the convention in the ensemble MCMC literature we call
those recipes of defining direction vectors, _moves_. Although the use of the
ensemble might seem equivalent to that of a sample covariance matrix in the
Adaptive Metropolis algorithm [77] the first has a higher information content
as it encodes both linear and non-linear correlations. Indeed, having an
ensemble of walkers allows for arbitrary many policies for choosing the
appropriate directions along which the walkers move in parameter space. As we
will shortly see, one of the choices (i.e. the Gaussian move, introduced later
in this Section) is indeed the slice sampling analogue of a covariance matrix.
However, other choices (i.e. Differential move or Global move) can take
advantage of the non-Gaussian nature of the ensemble distribution and thus
propose more informative moves. As will be discussed later in this section,
those advanced moves make no assumption of Gaussianity for the target
distribution. Furthermore, as we will show in the last part of this section,
the ensemble can also be easily parallelised.

Algorithm 16 Function to return a differential move direction vector.

1: function DifferentialMove($k$, $\mu$, $S$)

2: Draw two walkers $\mathbf{X_{l}}$, and $\mathbf{X_{m}}$ uniformly and
without replacement from the complementary ensemble $S$,

3: Compute direction vector $\bm{\eta}_{k}$ using Equation 10.6,

4: return $\bm{\eta}_{k}$

###### 45.2.2 Affine transformations and invariance

Affine invariance is a property of certain MCMC samplers first introduced in
the MCMC literature by [72]. An MCMC algorithm is said to be affine invariant
if its performance is invariant under the bijective mapping
$g:\mathbb{R}^{D}\rightarrow\mathbb{R}^{D}$ of the form
$\mathbf{Y}=A\mathbf{X}+b$ where $A\in\mathbb{R}^{D\times D}$ is a matrix and
$b\in\mathbb{R}^{D}$ is a vector. Linear transformations of this form are
called affine transformations and describe rotations, rescaling along specific
axes as well as translations in parameter space. Assuming that $\mathbf{X}$
has the probability density $p(\mathbf{X})$, then $\mathbf{Y}=A\mathbf{X}+b$
has the probability density

$p_{A,b}(\mathbf{Y})=p(A\mathbf{X}+b)\propto p(\mathbf{X})\,.$ (10.3)

Given a density $p$ as well as an MCMC transition operator $\mathcal{T}$ such
that $\mathbf{X}(t+1)=\mathcal{T}\big{(}\mathbf{X}(t);p\big{)}$ for any
iteration $t$ we call the operator $\mathcal{T}$ affine invariant if

$\mathcal{T}\big{(}A\mathbf{X}+b;p_{A,b}\big{)}=A\>\mathcal{T}\big{(}\mathbf{X};p\big{)}+b$
(10.4)

for $\forall A\in\mathbb{R}^{D\times D}$ and $\forall b\in\mathbb{R}^{D}$. In
case of an ensemble of walkers we define an affine transformation from
$\mathbb{R}^{ND}$ to $\mathbb{R}^{ND}$ as

$S=\\{\mathbf{X_{1}},\dots,\mathbf{X_{N}}\\}\xrightarrow{A,b}\\{A\mathbf{X_{1}}+b,\dots,A\mathbf{X_{N}}+b\\}\,.$
(10.5)

The property of affine invariance is of paramount importance for the
development of efficient MCMC methods. As we have discussed already, proposing
samples more frequently along certain directions can accelerate sampling by
moving further away in parameter space. Given that most realistic applications
are highly skewed or anisotropic and are characterised by some degree of
correlation between their parameters, affine invariant methods are an obvious
choice of a tool that can be used in order to achieve high levels of
efficiency.

###### 45.2.3 Differential move

The differential direction choice works by moving the walker $\mathbf{X}_{k}$
based on two randomly chosen walkers $\mathbf{X}_{l}$ and $\mathbf{X}_{m}$ of
the complementary ensemble $S_{[k]}=\\{\mathbf{X_{j}},\>\forall j\neq k\\}$
[69], see Figure 10.2 for a graphical explanation. In particular, we move the
walker $\mathbf{X}_{k}$ by slice sampling along the vector $\bm{\eta}_{k}$
defined by the difference between the walkers $\mathbf{X}_{l}$ and
$\mathbf{X}_{m}$. It is important to notice here that the vector
$\bm{\eta}_{k}$ is not a unit vector and thus carries information about both
the length scale and the optimal direction of movement. It will also prove to
be more intuitive to include the initial length scale $\mu$ in the definition
of the direction vector in the following way:

$\bm{\eta}_{k}=\mu\big{(}\mathbf{X}_{l}-\mathbf{X}_{m}\big{)}\,.$ (10.6)

The pseudocode for a function that, given the value of $\mu$ and the
complementary ensemble $S$, returns a differential direction vector
$\bm{\eta}_{k}$ is shown in Algorithm 16. Furthermore, the Differential move
is clearly affine invariant. Assuming that the distribution of the ensemble of
walkers follows the target distribution and the latter is highly elongated or
stretched along a certain direction then the proposed direction given by
equation 10.6 will share the same directional asymmetry.

Figure 10.2: The plot shows the differential direction move. Two walkers (red)
are uniformly sampled from the complementary ensemble (blue). Their positions
define the direction vector (solid black). The selected walker (magenta) then
moves by Slice Sampling along the parallel direction (dashed black).

###### 45.2.4 Gaussian move

The direction vector $\bm{\eta}_{k}$ can also be drawn from a normal
distribution with the zero mean and the covariance matrix equal to the sample
covariance of the complementary ensemble $S_{[k]}$,

$\mathbf{C}_{S}=\frac{1}{|S|}\sum_{j\in
S}\big{(}\mathbf{X}_{j}-\bar{\mathbf{X}}_{S}\big{)}\big{(}\mathbf{X}_{j}-\bar{\mathbf{X}}_{S}\big{)}^{t}\,.$
(10.7)

We chose to include the initial length scale $\mu$ in this definition as well:

$\frac{\bm{\eta}_{k}}{2\mu}\sim\mathcal{N}\big{(}\mathbf{0},\mathbf{C}_{S}\big{)}\,.$
(10.8)

The factor of $2$ is used so that the magnitude of the direction vectors are
consistent with those sampled using the differential direction choice in the
case of Gaussian-distributed walkers.

The pseudocode for a function that, given the value of $\mu$ and the
complementary ensemble $S$, returns a Gaussian direction vector
$\bm{\eta}_{k}$ is shown in Algorithm 17. See Figure 10.3 for a graphical
explanation of the method. Moreover, just like the Differential move, the
Gaussian move is also affine invariant. In the limit in which the number of
walkers is very large and the target distribution is normal, the first reduces
to the second. Alternatively, assuming that the distribution of walkers
follows the target distribution then the covariance matrix of the ensemble
would be the same as that of independently drawn samples from the target
density. Therefore any anisotropy characterising the target density would also
be present in the distribution of proposed directions given by equation 10.8.

Algorithm 17 Function to return a Gaussian Move direction vector.

1: function GaussianMove($k$, $\mu$, $S$)

2: Estimate sample covariance $\mathbf{C}_{S}$ of the walkers in the
complementary ensemble $S$ using Equation 10.7,

3: Sample
$\bm{\eta}_{k}/(2\mu)\sim\mathcal{N}\big{(}\mathbf{0},\mathbf{C}_{S}\big{)}$,

4: return $\bm{\eta}_{k}$

Figure 10.3: The plot shows the Gaussian direction move. A direction vector
(solid black) is sampled from the Gaussian-approximated distribution of the
walkers of the complementary ensemble (green). The selected walker (magenta)
then moves by Slice Sampling along the parallel direction (dashed black).

###### 45.2.5 Global move

ESS and its variations described so far (i.e. differential move, Gaussian
move) have as much difficulty traversing the low probability regions between
modes/peaks in multimodal distributions as most local MCMC methods (e.g.
Metropolis, Hamiltonian Monte Carlo, Slice Sampling, etc.). Indeed, multimodal
distributions are often the most challenging cases to sample from.
Fortunately, Ensemble Slice Sampling’s flexibility allows to construct
advanced moves which are specifically designed to handle multimodal cases even
in moderate to high dimensional parameter spaces. The _global move_ is such an
example.

We first fit a _Gaussian Mixture_ to the distribution of the walkers of the
complementary ensemble $S_{[k]}$ using _Variational Inference_. To avoid
defining the number of components of the Gaussian Mixture we use a _Dirichlet
process_ as the prior distribution for the Gaussian Mixture weights111To this
end we use the Scikit-Learn implementation of the Dirichlet process Gaussian
mixture. [73]. The exact details of the construction of the Dirchlet process
Gaussian mixture (DPGM) are beyond the scope of this work and we direct the
reader to [73] and [18] for more details. One of the major benefits of fitting
the DPGM using variational inference compared to the expectation–maximisation
(EM) algorithm [33] that is often used is the improved stability. In
particular, the use of priors in the variational Bayesian treatment guarantees
that Gaussian components do not collapse into specific data points. This
regularisation due to the priors leads to component covariance matrices that
do not diverge even when the number of data points (i.e. walkers in our case)
in a component is lower than the number of dimensions. In our case, this means
that even if the number of walkers located in a mode of the target
distribution is small DPGM would still identify that mode correctly. In such
cases, the covariance of the component that corresponds to that mode would be
over–estimated. This however does not affect the performance of the Global
move as the latter does not rely on exact estimates of the component
covariance matrices.222Indeed the covariance matrix of a component only enters
through equation 10.10 but then it is re–scaled by the factor $\gamma$.

In practice, we recommend using more than the minimum number of walkers in
cases of multimodal distributions (e.g. at least two times as many in bimodal
cases). We found that the computational overhead introduced by the variational
fitting of the DPGM is negligible compared to the computational cost of the
evaluation of the model and posterior distribution in common problems in
physics, astrophysics and cosmology. Indeed the cost is comparable, and only a
few times higher than the Differential or Gaussian move. The reason for that
is the relatively small number of walkers (i.e. $\mathcal{O}(10-10^{3})$) that
simplifies the fitting procedure.

Once fitting is done, we have a list of the means and covariance matrices of
the components of the Gaussian Mixture. As the ensemble of walkers traces the
structure of the target distribution, we can use the knowledge of the means
and covariance matrices of the Gaussian Mixture to construct efficient
direction vectors. Ideally, we prefer direction vectors that connect different
modes. This way, the walkers will be encouraged to move along those directions
that would otherwise be very unlikely to be chosen.

We uniformly select two walkers of the complementary ensemble and identify the
Gaussian components to which they belong, say $i$ and $j$. There are two
distinct cases and we will treat them as such. In case A, $i=j$, meaning that
the selected walkers originate from the same component. In case B, $i\neq j$,
meaning that the two walkers belong to different components and thus probably
different peaks of the target distribution.

As we will show next, only in case B, we can define a direction vector that
favors mode-jumping behaviour. In case A, we can sample a direction vector
from the Gaussian component that the two select walkers belong to333In
practice we use uniformly sample two walkers from the list of walkers that
DPGM identified in that mode. This step removes any dependency on covariance
matrix estimates.:

$\frac{\bm{\eta}_{k}}{2\mu}\sim\mathcal{N}\big{(}\bm{0},\bm{C}_{i=j}\big{)}\,,$
(10.9)

where $\bm{C}_{i=j}$ is the covariance matrix of the ith (or equivalently jth)
component. Just as in the Gaussian move, the mean of the proposal distribution
is zero so that we can interpret $\bm{\eta}$ as a direction vector.

In case B, where the two selected walkers belong to different components,
$i\neq j$, we will follow a different procedure to facilitate long jumps in
parameter space. We will sample two vectors, one from each component:

$\bm{\eta}_{k,n}\sim\mathcal{N}\big{(}\bm{\mu}_{n},\gamma\bm{C}_{n}\big{)}\,,$
(10.10)

for $n=i$ or $n=j$. Here, $\bm{\mu}_{n}$ is the mean of the nth component and
$\bm{C}_{n}$ is its covariance matrix. In practice, we also re-scale the
covariance by a factor of $\gamma=0.001$, which results in direction vectors
with lower variance in their orientation. $\gamma<1$ ensures that the chosen
direction vector is close to the vector connecting the two peaks of the
distribution. Finally, the direction vector will be defined as:

$\bm{\eta}_{k}=2\big{(}\bm{\eta}_{k,i}-\bm{\eta}_{k,j}\big{)}\,.$ (10.11)

The factor of $2$ here is chosen to better facilitate mode-jumping. There is
also no factor of $\mu$ in the aforementioned expression since in this case
there is no need for the scale factor to be tuned.

The pseudocode for a function that, given the complementary ensemble $S$,
returns a Global direction vector $\bm{\eta}_{k}$ is shown in Algorithm 18.
See Figure 10.4 for a graphical explanation of the method. It should be noted
that for the global move to work at least one walker needs to be present on
each well separated mode.

Algorithm 18 Function to return a global move direction vector.

1: function GlobalMove($k$, $\mu$, $S$)

2: Fit Dirichlet process Gaussian mixture (DPGM) to the complementary ensemble
$S_{[k]}$,

3: If $N$ is the number of components of the DPGM then select two components
$i,j$ uniformly such that $i\neq j$,

4: if $i=j$ then

5: Sample
$\bm{\eta}_{k}/(2\mu)\sim\mathcal{N}\big{(}\bm{0},\bm{C}_{i=j}\big{)}$,

6: else

7: Sample
$\bm{\eta}_{k,n}\sim\mathcal{N}\big{(}\bm{\mu}_{n},\gamma\bm{C}_{n}\big{)}$
for $n=i,j$,

8: Compute direction vector $\bm{\eta}_{k}$ using Equation 10.11,

9: end if

10: return $\bm{\eta}_{k}$

Figure 10.4: The plot shows the global direction move assuming that the
uniformly selected pair of walkers of the complementary ensemble belongs to
different components (blue and green). A position (red) is sampled from each
component (using the re-scaled by $\gamma$ covariance matrix). Those two
points (red) define the direction vector (black) connecting the two modes
(blue and green). The selected walker (magenta) then moves by slice sampling
along the parallel direction (dashed).

Here we introduced three general and distinct moves that can be used in a
broad range of cases. In general, the global move requires a higher number of
walkers than the differential or Gaussian move in order to perform well. We
found that the differential and Gaussian moves are good choices for most
target distributions whereas the global move is only necessary in highly
dimensional and multimodal cases. One can use the information in the
complementary ensemble to construct more moves tailor-made for specific
problems. Such additional moves might include Kernel Density Estimation or
Clustering methods and as long as the information used comes from the
complementary ensemble (and not from the walker that would be updated) the
detailed balance is preserved.

###### 45.2.6 Parallelizing the ensemble

Instead of evolving the ensemble by moving each walker in turn we can do this
in parallel. A naive implementation of this would result in a subtle violation
of detailed balance. We can avoid this by splitting the ensemble into two sets
of walkers [50] of $n_{\text{Walkers}}/2$ each. We can now update the
positions of all the walkers in the one set in parallel along directions
defined by the walkers of the other set (the complementary ensemble). Then we
can perform the same procedure for the other set. In accordance with equation
10.2, the stationary distribution of the split ensemble would be

$P(\mathbf{X_{1}},\dots,\mathbf{X_{N}})=\prod_{k=1}^{N/2}p(\mathbf{X_{k}})\prod_{k=1+N/2}^{N}p(\mathbf{X_{k}})\,.$
(10.12)

The method generates samples from the target distribution by simulating a
Markov chain which leaves this product distribution invariant. The transition
operator $\mathcal{T}_{1}$ that updates the walkers of the first set (i.e.
$k=1,\dots,N/2$) uses the walkers of the complementary ensemble (i.e.
$k=1+N/2,\dots,N$) and vice versa for the transition operator
$\mathcal{T}_{2}$ that acts on the second set. In the context of ESS the
aforementioned transition operators correspond to a single iteration of
Algorithm 19 coupled with one of the moves (e.g. Differential move).

It follows from the ensemble splitting technique that the maximum number of
CPUs used without any of them being idle is equal to the total number of
walkers updated concurrently, that is $n_{\text{Walkers}}/2$. We will also
verify this empirically in Section 46. Of course, this does not mean that if
there are more CPUs available they cannot be used as we can always increase
the size of the ensemble to match the available CPUs.

Combining this technique with the stochastic approximation solution of
Subsection 45.1 and the choices (moves) of direction and ensemble-splitting
technique of this subsection leads to the Ensemble Slice Sampling method of
Algorithm 19444Perhaps a small detail, but we have included the length scale
in the definition of the direction vector $\mathbf{\eta}$ and therefore it
does not appear in the definition of the $(L,R)$ interval.. Of course, another
move (e.g. Gaussian, global) can be used instead of the differential move in
Algorithm 19. Finally, the minimum number of walkers used should be twice the
number of parameters. Using fewer walkers than that could lead to erroneous
sampling from a lower dimensional parameter space [169].

In general, parallelizing a slice sampler is not trivial (e.g. as it is for
Metropolis) because each update requires an unknown number of probability
density evaluations. However, because of the affine invariance (i.e.
performance unaffected by linear correlations) induced by the existence of the
ensemble, all iterations require on average the same number of probability
density evaluations (i.e. usually $5$ if the stochastic approximation for the
length scale $\mu$ is used). Therefore, the parallelization of Ensemble Slice
Sampling is very effective in practice. Furthermore, the benefit of having
parallel walkers instead of parallel independent chains (e.g. such as in
Metropolis sampling) is clear, the walkers share information about the
covariance structure of the distribution thus accelerating mixing.

Algorithm 19 Single Iteration $t$ of Ensemble Slice Sampling.

1: Given $t$, $f$, $\mu^{(t)}$, $S_{[0]}$, $S_{[1]}$, $M^{\rm adapt}$:

2: Initialise $N_{e}^{(t)}=0$ and $N_{c}^{(t)}=0$,

3: for $i=0,1$ do

4: for $k=1,...,N/2$ do

5: $k\leftarrow k+iN/2$

6: Compute direction vector $\bm{\eta}_{k}\leftarrow$ DifferentialMove($k$,
$\mu^{(t)}$, $S_{[i]}$)

7: Sample $Y\sim\text{Uniform}(0,f(\mathbf{X_{k}}^{(t)}))$

8: Sample $U\sim\text{Uniform}(0,1)$

9: Set $L\leftarrow-U$, and $R\leftarrow L+1$

10: while $Y<f(L)$ do

11: $L\leftarrow L-1$

12: $N_{e}^{(t)}\leftarrow N_{e}^{(t)}+1$

13: end while

14: while $Y<f(R)$ do

15: $R\leftarrow R+1$

16: $N_{e}^{(t)}\leftarrow N_{e}^{(t)}+1$

17: end while

18: while True do

19: Sample $X^{\prime}\sim\text{Uniform}(L,R)$

20: Set $Y^{\prime}\leftarrow f(X^{\prime}\bm{\eta}_{k}+\mathbf{X_{k}}^{(t)})$

21: if $Y<Y^{\prime}$ then

22: break

23: end if

24: if $X^{\prime}<0$ then

25: $L\leftarrow X^{\prime}$

26: $N_{c}^{(t)}\leftarrow N_{c}^{(t)}+1$

27: else

28: $R\leftarrow X^{\prime}$

29: $N_{c}^{(t)}\leftarrow N_{c}^{(t)}+1$

30: end if

31: end while

32: Set $\mathbf{X_{k}}^{(t+1)}\leftarrow
X^{\prime}\bm{\eta}_{k}+\mathbf{X_{k}}^{(t)}$

33: end for

34: end for

35: $\mu^{(t+1)}\leftarrow$ TuneLengthScale($t$, $\mu^{(t)}$, $N_{e}^{(t)}$,
$N_{c}^{(t)}$, $M^{adapt}$),

#### 46 Empirical evaluation

To empirically evaluate the sampling performance of the Ensemble Slice
Sampling algorithm we perform a series of tests. In particular, we compare its
ability to sample from two demanding target distributions, namely the
_autoregressive process of order 1_ and the _correlated funnel_ , against the
Metropolis and Standard Slice Sampling algorithms. The Metropolis’ proposal
scale was tuned to achieve the optimal acceptance rate, whereas the initial
length scale of Standard Slice Sampling was tuned using the stochastic scheme
of Algorithm 15. Ensemble Slice Sampling significantly outperforms both of
them. These tests help establish the characteristics and advantages of
Ensemble Slice Sampling. Since our objective was to develop a gradient-free
black-box method we then proceed to compare Ensemble Slice Sampling with a
list of gradient-free ensemble methods such as _Affine Invariant Ensemble
Sampling_ (AIES), _Differential Evolution Markov Chain_ (DEMC) and _Kernel
Density Estimate Metropolis_ (KM) on a variety of challenging target
distributions. Moreover, we are also interested in assessing the convergence
rate of the length scale $\mu$ during the first iterations as well as the
parallel scaling of the method in the presence of multiple CPUs. Unless
otherwise specified we use the differential move for the tests. Unlike ESS
that has an acceptance rate of $1$, AIES’s and DEMC’s acceptance rate is
related to the number of walkers. For that reason, and for the sake of a fair
comparison, we made sure the selected number of walkers in all examples would
yield the optimal acceptance rate for AIES and DEMC. As we will discuss
further in Section 47 it makes sense to increase the number of walkers in
cases of multimodal distributions or strong non-linear correlations. In
general though, we recommend using the minimum number of walkers (i.e. twice
the number of dimensions) as the default choice and increase it only if it is
required by a specific application. For more rules and heuristics about the
initialisation and number of walkers we direct the interested reader to
Section 47.

##### 46.1 Performance tests

###### 46.1.1 Autoregressive process of order 1

In order to investigate the performance of ESS. in high dimensional and
correlated scenarios we chose a highly correlated Gaussian as the target
distribution. More specifically, the target density is a discrete-time
_autoregressive process of order 1_ , also known as AR(1). This particular
target density is ideally suited for benchmarking MCMC algorithms since the
posterior density in many scientific studies often approximates a correlated
Gaussian. Apart from that, the AR(1) is commonly used as a prior for time-
series analysis.

The AR(1) distribution of a random vector $\bm{X}=(X_{1},...,X_{N})$ is
defined recursively as follows:

$\begin{split}X_{1}\sim&\;\mathcal{N}(0,1)\,,\\\
X_{2}|X_{1}\sim&\;\mathcal{N}(\alpha X_{1},\beta^{2})\,,\\\ &\vdots\\\
X_{N}|X_{N-1}\sim&\;\mathcal{N}(\alpha X_{N-1},\beta^{2})\,,\end{split}$
(10.13)

where the parameter $\alpha$ controls the degree of correlation between
parameters and we chose it to be $\alpha=0.95$. We set
$\beta=\sqrt{1-\alpha^{2}}$ so that the marginal distribution of all
parameters is $\mathcal{N}(0,1)$. We also set the number of dimensions to
$N=50$.

Figure 10.5: The plots compare the 1-sigma and 2-sigma contours generated by
the optimised random-walk Metropolis (left), Standard Slice (centre) and
Ensemble Slice Sampling (right) methods to those obtained by Independent
Sampling (blue) for the AR(1) distribution. All samplers used the same number
of probability density evaluations, $3\times 10^{5}$. Only the first two
dimensions are shown here. Table 2: The table shows a comparison of the
optimally tuned Metropolis, Standard Slice, and Ensemble Slice Sampling with
the differential move (ESS-D) and the Gaussian move (ESS-G) respectively in
terms of the integrated autocorrelation time (IAT) and the number of effective
samples per evaluation of the probability density (efficiency) multiplied by
$10^{4}$. These metrics are formally defined in Appendix 49. The target
distributions are the 50–dimensional autoregressive process of order 1 and the
25–dimensional correlated funnel distribution. The total number of iterations
was set to $10^{7}$.

[0.75pt] | Metropolis | Slice | ESS-D | ESS-G
---|---|---|---|---
[0.5pt] Autoregressive process of order 1 |
[0.5pt] IAT | 4341 | 2075 | $\mathbf{111}$ | $\mathbf{107}$
efficiency | 2.3 | 1.0 | $\mathbf{17.5}$ | $\mathbf{17.8}$
[0.5pt] Correlated funnel distribution |
[0.5pt] IAT | - | 3905 | $\mathbf{129}$ | $\mathbf{141}$
efficiency | - | 0.5 | $\mathbf{15.3}$ | $\mathbf{14.0}$
[0.75pt] |  |  |  |

For each method, we measured the mean _integrated autocorrelation time_ (IAT),
and the number of effective samples per evaluation of the probability density
function, also termed _efficiency_ (see Appendix 49 for details). For this
test we ran the samplers for $10^{7}$ iterations. In this example we used the
minimum number of walkers (i.e. 100 walkers) for ESS and the equivalent number
of probability evaluations for Metropolis and Slice Sampling with each walker
initialised at a position sampled from the distribution $\mathcal{N}(0,1)$.
The results are presented in Table 46.1.1. The chain produced by Ensemble
Slice Sampling has a significantly shorter IAT ($20-40$ times) compared to
either of the other two methods. Furthermore, Ensemble Slice Sampling, with
either Differential or Gaussian move, generates an order of magnitude greater
number of independent samples per evaluation of the probability density. In
this example the Differential and Gaussian moves have achieved almost
identical IAT values and efficiencies.

To assess the mixing rate of Ensemble Slice Sampling, we set the maximum
number of probability density evaluations to $3\times 10^{5}$ and show the
results in Figure 10.5. We compare the results of Ensemble Slice Sampling with
those obtained via the optimally tuned Metropolis and Standard Slice Sampling
methods. Ensemble Slice Sampling significantly outperforms both of them, being
the only one with a chain resembling the target distribution in the chosen
number of probability evaluations.

###### 46.1.2 Correlated funnel

The second test involves a more challenging distribution, namely the
correlated funnel distribution adapted from [133]. The funnel, tornado like,
structure is common in Bayesian hierarchical models and possesses
characteristics that render it a particularly difficult case. The main
difficulty originates from the fact that there is a region of the parameter
space where the volume of the region is low but the probability density is
high, and another region where the opposite holds.

Suppose we want to sample an N–dimensional vector $\bm{X}=(X_{1},...,X_{N})$
from the correlated funnel distribution. The marginal distribution of $X_{1}$
is Gaussian with mean zero and unit variance. Conditional on a value of
$X_{1}$, the vector $\bm{X}_{2-N}=(X_{2},...,X_{N})$ is drawn from a Gaussian
with mean zero and a covariance matrix in which the diagonal elements are
$\exp(X_{1})$, and the non-diagonal equal to $\gamma\exp(X_{1})$. If
$\gamma=0$, the parameters $X_{2}$ to $X_{N}$ conditional on $X_{1}$ are
independent and the funnel distribution resembles the one proposed by [133].
The value of $\gamma$ controls the degree of correlation between those
parameters. When $\gamma=0$ the parameters are uncorrelated. For the following
test we chose this to be $\gamma=0.95$. We set the number of parameters $N$ to
$25$.

Using $10^{7}$ iterations, we estimated the IAT and the efficiency of the
algorithms for this distribution as shown in Table 46.1.1. Just like in the
AR(1) case we used the minimum number (i.e. 50) of walkers for ESS with each
walker initialised at a position sampled from the distribution
$\mathcal{N}(0,1)$. Since the optimally-tuned Metropolis fails to sample from
this particular distribution, we do not quote any results. The Metropolis
sampler is unable to successfully explore the region of parameter space with
negative $X_{1}$ values. The presence of strong correlations renders the
Ensemble Slice Sampler $30$ times more efficient than the Standard Slice
Sampling algorithm on this particular example. In this example, the
Differential move outperforms the Gaussian move in terms of efficiency, albeit
by a small margin. In general, we expect the former to be more flexible than
the latter since it makes no assumption about the Gaussianity of the target-
distribution and recommend it as the default configuration of the algorithm.

Figure 10.6: The plots compare the 1-sigma and 2-sigma contours generated by
the optimised random-walk Metropolis (left), Standard Slice (centre) and
Ensemble Slice Sampling (right) methods to those obtained by Independent
Sampling (blue) for the correlated funnel distribution. All samplers used the
same number of probability density evaluations, $3\times 10^{5}$. Only the
first two dimensions are shown here.

To assess the mixing rate of the algorithm on this demanding case, we set the
maximum number of evaluations of the probability density function to $3\times
10^{5}$. As shown in Figure 10.6, the Ensemble Slice Sampling is the only
algorithm out of the three whose outcome closely resembles the target
distribution. The results of Metropolis were incorrect for both, the limited
run with $3\times 10^{5}$ iterations and the long run with $10^{7}$
iterations. In particular, the chain produced using the Metropolis method
resemble a converged chain but in fact it is biased in favour of positive
values of $x_{1}$. The problem arises because of the vanishing low probability
of accepting a point with highly negative value of $x_{1}$. This indicates the
inability of Metropolis to handle this challenging case. For a more detailed
discussion of this problem we direct the reader to Section 8 of [133]. In
general, the correlated funnel is a clear example of a distribution in which a
single Metropolis proposal scale is not sufficient for all the sampled regions
of parameter space. The locally adaptive nature of ESS solves this issue.

##### 46.2 Comparison to other ensemble methods

So far we have demonstrated Ensemble Slice Sampling’s performance in simple,
yet challenging, target distributions. The tests performed so far demonstrate
ESS’s capacity to sample efficiently from highly correlated distributions
compared with standard methods such as Metropolis and Slice Sampling. Although
the use of Metropolis and Slice Sampling is common, these methods are not
considered to be state-of-the-art. For this reason, we will now compare ESS
with state-of-the-art gradient-free ensemble MCMC methods.

By far, the two most popular choices555For instance, in the fields of
Astrophysics and Cosmology where most models are not differentiable and
gradient methods (e.g. Hamiltonian Monte Carlo or NUTS) are not applicable the
default choice is the Affine-Invariant Ensemble Sampler (AIES) [72] as
implemented in emcee. of gradient-free ensemble methods are the Affine-
Invariant Ensemble Sampling (AIES) [72] method and the Differential Evolution
Monte Carlo (DEMC) [169] algorithm supplemented with a Snooker update [170].

In cases of strongly multimodal target distributions we will also test our
method against Sequential Monte Carlo666As there are many different flavours
of SMC, we decided to use the one implemented in PyMC3 which utilises
importance sampling, simulated annealing and Metropolis sampling. (SMC) [115,
32] and Kernel Density Estimate Metropolis (KM) [43] which are particle
methods specifically designed to handle strongly multimodal densities.

###### 46.2.1 Ring distribution

Although, all three of the compared methods (i.e. ESS, AIES, DEMC) are affine
invariant and thus unaffected by linear correlations, they do however differ
significantly in the way they handle non-linear correlations. In particular,
only Ensemble Slice Sampling (ESS) is locally adaptive because of its
stepping-out procedure and therefore able to handle non-linear correlations
efficiently.

To illustrate ESS’s performance in a case of strong non-linear correlations we
will use the 16–dimensional ring distribution defined by:

$\begin{split}\ln\mathcal{L}=&-\Bigg{[}\frac{(x_{n}^{2}+x_{1}^{2}-a)^{2}}{b}\Bigg{]}^{2}\\\
&-\sum_{i=1}^{n-1}\Bigg{[}\frac{(x_{i}^{2}+x_{i+1}^{2}-a)^{2}}{b}\Bigg{]}^{2}\,,\end{split}$
(10.14)

where $a=2$, $b=1$ and $n=16$ is the total number of parameters. We also set
the number of walkers to be $64$ and run the samplers for $10^{7}$ steps
discarding the first half of the chains. Here we followed the heuristics
discussed at the beginning of this section and increased the number of walkers
from the minimum of $2\times 16$ to $4\times 16$ due to the presence of strong
non-linear correlations in order to achieve the optimal acceptance rate for
AIES and DEMC. The number of iterations is large enough for all samplers to
converge and provide accurate estimates of the autocorrelation time.

The results are shown in Table LABEL:tab_ess:table2 and verify that ESS’
performance is an order of magnitude better than that of the other methods.

Table 3: The table shows a comparison of the Affine Invariant Ensemble
Sampling (AIES), Differential Evolution Markov Chain (DEMC), and Ensemble
Slice Sampling methods in terms of the integrated autocorrelation time (IAT)
and the number of effective samples per evaluation of the probability density
(efficiency) multiplied by $10^{5}$. These metrics are formally defined in
Appendix 49. The target distributions are the 16–dimensional ring
distribution, the 10–dimensional Gaussian shells distribution and the
13–dimensional hierarchical Gaussian process regression distribution. In all
cases the total number of iterations was set to $10^{7}$. It should be noted
that in the case of the Gaussian shells the global move was used instead of
the differential move.

###### 46.2.3 Hierarchical Gaussian process regression

To illustrate ESS’s performance in a real-world example we will use a
modelling problem concerning the concentration of $CO_{2}$ in the atmosphere
adapted from Chapter 5 of [146]. The data consist of monthly measurements of
the mean $CO_{2}$ concentration in the atmosphere measured at the _Mauna Loa
Observatory_ [100] in _Hawaii_ since 1958. Our goal is to model the
concentration of $CO_{2}$ as a function of time. To this end, we will employ a
_hierarchical Gaussian process_ model with a composite covariance function
designed to take care of the properties of the data. In particular, the
covariance function (kernel) is the sum of following four distinct terms:

$k_{1}(r)=\theta_{1}^{2}\exp\left(-\frac{r^{2}}{2\theta_{2}}\right)\,,$
(10.17)

where $r=x-x^{\prime}$ that describes the smooth trend of the data,

$k_{2}(r)=\theta_{3}^{2}\exp\left[-\frac{r^{2}}{2\theta_{4}}-\theta_{5}\sin^{2}\left(\frac{\pi
r}{\theta_{6}}\right)\right]\,,$ (10.18)

that describes the seasonal component,

$k_{3}(r)=\theta_{7}^{2}\left[1+\frac{r^{2}}{2\theta_{8}\theta_{9}}\right]^{-\theta_{8}}\,,$
(10.19)

which encodes medium-term irregularities, and finally:

$k_{4}(r)=\theta_{10}^{2}\exp\left(-\frac{r^{2}}{2\theta_{11}}\right)+\theta_{12}^{2}\delta_{ij}\,,$
(10.20)

that describes the noise. We also fit the mean of the data, having in total 13
parameters to sample.

We sample this target distribution using $36$ walkers for $10^{7}$ iterations
and we discard the first half of the chains. The number of walkers that was
used corresponds to $1.5$ times the minimum number. We found that this value
results in the optimal acceptance rate for AIES and DEMC. For this example we
use the differential move of ESS. The results are presented in Table
LABEL:tab_ess:table2. The integrated autocorrelation time of ESS is $2$ orders
of magnitude lower than that of the other methods and its efficiency is more
than an order of magnitude higher. The performance is weakly sensitive to the
choice of the number of walkers.

###### 46.2.4 Bayesian object detection

Another real world example with many applications in the field of _astronomy_
is _Bayesian object detection_. The following model adapted from [45] can be
used with a few adjustments to detect astronomical objects in telescope images
often hidden in background noise.

We assume that the 2D circular objects present in the image are described by
the Gaussian profile:

$\mathbf{G}(x,y;\bm{\theta})=A\exp\bigg{[}-\frac{(x-X)^{2}+(y-Y)^{2}}{2R^{2}}\bigg{]}\,,$
(10.21)

where $\mathbf{\theta}=(X,Y,A,R)$ are parameters that define the coordinate
position, the amplitude and the size of the object, respectively. Then the
data can be described as:

$\mathbf{D}=\mathbf{N}+\sum_{i=1}^{n_{\text{Obj}}}\mathbf{G}(\bm{\theta_{i}})\,,$
(10.22)

where $n_{\text{Obj}}$ is the number of objects in the image and $\mathbf{N}$
is an additive Gaussian noise term.

Assuming a $200\times 200$ pixel-wide image, we can create a simulated dataset
by sampling the coordinate positions $(X,Y)$ of the objects from
$\mathcal{U}(0,200)$ and their amplitude $A$ and size $R$ from
$\mathcal{U}(1,2)$ and $\mathcal{U}(3,7)$, respectively. We sample
$n_{\text{Obj}}=8$ objects in total. Finally, we sample the noise $\mathbf{N}$
from $\mathcal{N}(0,4)$. In practice we create a dataset of $100$ such images
and one such example is shown in Figure 10.7. Notice that the objects are
hardly visible as they are obscured by the background noise, this makes the
task of identifying those objects very challenging.

Figure 10.7: The plot shows a simulated image used in the Bayesian object
detection exercise. There are $8$ circular objects included here. As the
objects are hardly visible due to the background noise their centres are
marked with red stars.

Following the construction of the simulated dataset, the posterior probability
density function is defined as:

$P(\bm{\theta}|\mathbf{D})\propto\exp\bigg{\\{}\frac{[\mathbf{G}(\bm{\theta})-\mathbf{D}]^{2}}{2\sigma^{2}}\bigg{\\}}P(\bm{\theta})\,,$
(10.23)

where $\sigma=2$ is the standard deviation of the $\mathbf{N}$ noise term. The
prior $P(\bm{\theta})$ can be decomposed as the product of prior distributions
of $X$, $Y$, $A$, and $R$. We used uniform priors for all of these parameters
with limits $(0,200)$ for $X$ and $Y$, $(1,2)$ for $A$, and $(2,9)$ for $R$.
It is important to mention here that the posterior does not include any prior
information about the exact or maximum number of objects in the data. In that
sense, the sampler is agnostic about the exact number, positions and
characteristics (i.e. amplitude and size) of the objects that it seeks to
detect.

We sampled the posterior distribution using $200$ walkers (initialised from
the prior distribution) for each image in our dataset (i.e. 100 images in
total) using Ensemble Slice Sampling (ESS), Affine Invariant Ensemble Sampling
(AIES), and Differential Evolution Markov Chain (DEMC). Although the posterior
distribution is multimodal (i.e. $8$ modes) we used the differential move
since the number of dimensions is low and there is no reason to use more
sophisticated moves like the global move. We used a large enough ensemble of
walkers due to the potential presence of multiple modes so that all three
samplers are able to resolve them.

We ran each sampler for $10^{4}$ iterations in total and we discarded the
first half of the chains. We found that, on average for the 100 images, ESS
identifies correctly $7$ out of $8$ objects in the image, whereas AIES and
DEMC identify $4$ and $5$, respectively.

In cases where the objects are well-separated ESS often identifies correctly
$8$ out of $8$. Its accuracy falls to $7/8$ in cases where two of the objects
are very close to each other or overlap. In those cases ESS identifies the
merged object as a single object. In this context, by identification of an
object, we mean that at least one walker has sampled the posterior mode which
corresponds to that object.

###### 46.2.5 Gaussian Mixture

One strengths of ESS is its ability to sample from strongly multimodal
distributions in high dimensions. To demonstrate this, we will utilise a
Gaussian Mixture of two components centred at $\mathbf{-0.5}$ and
$\mathbf{+0.5}$ with standard deviation of $\mathbf{0.1}$. We also put $1/3$
of the probability mass in one mode and $2/3$ in the other.

We first set this distribution at $10$ dimensions and we sample this using
$80$ walkers for $10^{5}$ steps. The distance between the two modes in this
case is approximately $32$ standard deviations. We then increase the number of
dimensions to $50$ and we sample it using $400$ walkers for $10^{5}$
iterations. In this case, the actual distance between the two modes is
approximately $71$ standard deviations. The total number of iterations was set
to $10^{7}$ for all methods but the SMC.

This problem consists of two, well separated, modes and thus requires using at
least twice the minimum number of walkers (i.e. at least 40 for the
10–dimensional case and 200 for the 50–dimensional one). Although the
aforementioned configuration was sufficient for ESS to provide accurate
estimates, we opted instead for twice that number (i.e. 80 walkers for the
10–dimensional cases and 400 for the 50–dimensional one) in order to satisfy
the requirements of the other samplers, mainly the Kernel Density Estimate
Metropolis (KM), but also AIES and DEMC. For the Sequential Monte Carlo (SMC)
sampler we used $2000$ and $20000$ independent chains for the low and high
dimensional case respectively. The temperature ladder that interpolates
between the prior and posterior distribution was chosen adaptively
guaranteeing an effective sample size of $90\%$ the physical size of the
ensemble. Our implementation of SMC was based on that of PyMC3 using an
independent Metropolis mutation kernel.

The results for the 10–dimensional and 50–dimensional cases are plotted in
Figures 10.8 and 10.9, respectively. In the 10–dimensional case, both ESS
(differential and global move) and SMC managed to sample from the target
whereas AIES, DEMC and KM failed to do so. In the 50–dimensional case, only
the Ensemble Slice Sampling with the global move manages to sample correctly
from this challenging target distribution. In practice $\text{ESS}_{G}$ is
able to handle similar cases in even higher number of dimensions and with more
than $2$ modes.

Figure 10.8: The plot compares the results of 6 samplers, namely Sequential
Monte Carlo (SMC, red), Affine-Invariant Ensemble Sampling (AIES, yellow),
Differential Evolution Markov Chain (DEMC, purple), Kernel Density Estimate
Metropolis (KM, orange), Ensemble Slice Sampling using the differential move
(ESS, green), and Ensemble Slice Sampling using the global move (ESS, blue).
The target distribution is a 10–dimensional Gaussian Mixture. The figure shows
the 1D marginal distribution for the first parameter of the 10. Figure 10.9:
The plot compares the results of 6 samplers, namely Sequential Monte Carlo
(SMC, red), Affine-Invariant Ensemble Sampling (AIES, yellow), Differential
Evolution Markov Chain (DEMC, purple), Kernel Density Estimate Metropolis (KM,
orange), Ensemble Slice Sampling using the differential move (ESS, green), and
Ensemble Slice Sampling using the global move (ESS, blue). The target
distribution is a 50–dimensional Gaussian Mixture.The figure shows the 1D
marginal distribution for the first parameter of the 50.

##### 46.3 Convergence of the Length Scale $\mu$

Figure 10.10 plots the convergence of the length scale during the first 20
iterations. The target distribution in this example is a 20–dimensional
correlated normal distribution. The length scale $\mu$ was initialised from a
wide range of possible values. Adaptation is significantly faster when the
initial length scale is larger than the optimal one rather than smaller.
Another benefit of using a larger initial estimate would be the reduced number
of probability evaluations during the first iterations. This is due to the
fact that the shrinking procedure is generally faster than the stepping-out
procedure.

Figure 10.10: The plot shows the adaptation of the length scale $\mu$ as a
function of the number of iterations and starting from a wide range of initial
values. Each trace is an independent run and the y-axis shows the value of
$\mu$ divided by the final value of $\mu$. The target distribution in this
example is a 20–dimensional correlated normal distribution. Starting from
larger $\mu$ values leads to significantly faster adaptation.

##### 46.4 Parallel Scaling

By construction, Ensemble Slice Sampling can be used in parallel computing
environments by parallelising the ensemble of walkers as discussed in Section
45.2. The maximum number of CPUs used without any of them being idle is equal
to the size of complementary ensemble, $n_{\text{Walkers}}/2$. In order to
verify this empirically and investigate the scaling of the method for any
number of CPUs, we sampled a 10–dimensional Normal distribution for $10^{5}$
iterations with varying number of walkers. The results are plotted in Figure
10.11. We sampled the aforementioned distribution multiple times in order to
get estimates of the confidence integrals shown in Figure 10.11. The required
time to do the pre-specified number of iterations scales as
$\mathcal{O}(1/n_{\text{CPUs}})$ as long as $n_{\text{CPUs}}\leq
n_{\text{Walkers}}/2$. This result does not depend on the specific
distribution. We can always use all the available CPUs by matching the size of
the complementary ensemble (i.e. half the number of walkers) to the number of
CPUs.

Figure 10.11: The plot shows the time $t_{f}$ required for ESS to complete a
pre-specified number of iterations as a function of the ratio of the number of
available CPUs $n_{\rm CPUs}$ to the total number of walkers $n_{\rm
Walkers}$. The results are normalised with respect to the single CPU case
$t_{1}$. The method scales as $\mathcal{O}(1/n_{\text{CPUs}})$ as long as
$n_{\text{CPUs}}\leq n_{\text{Walkers}}/2$ (dashed line). The shaded areas
show the $2-\sigma$ intervals.

#### 47 Discussion

In Section 46 we provided a quantitative comparison of the efficiency of
Ensemble Slice Sampling compared to other methods. In this Section we will
provide some qualitative arguments to informally demonstrate the advantages of
Ensemble Slice Sampling over other methods. Furthermore, we will briefly
discuss some general aspects of the algorithm and place our work in the
context of other related algorithms.

After the brief adaptation period is over and the length scale $\mu$ is fixed,
the Ensemble Slice Sampling algorithm performs on average $5$ evaluations of
the probability density per walker per iteration, assuming that either the
differential or Gaussian move is used. This is in stark contrast with
Metropolis-based MCMC methods that perform $1$ evaluation of the probability
density per iteration. However, the non-rejection nature of Ensemble Slice
Sampling more than compensates for the higher number of evaluations as shown
in Section 46, thus yielding a very efficient scheme.

One could think of the number of walkers as the only free hyperparameter of
Ensemble Slice Sampling. However, choosing the number of walkers is usually
trivial. As we mentioned briefly at the end of Section 45, there is a minimum
limit to that number. In particular, in order for the method to be ergodic,
the ensemble should be made of at least $2\times D$ walkers777The reason that
the minimum limit is $2\times D$ instead of $D+1$ has to do with the ensemble
splitting procedure that we introduced in order to make the method parallel.
Splitting the ensemble into two equal parts means that each walker is updated
based on the relative displacements of half the ensemble., where $D$ is the
number of dimensions of the problem. Assuming that the initial relative
displacements of the walkers span the parameter space (i.e. they do not belong
to a lower-than-$D$-dimensional space) the resulting algorithm would be
ergodic. As shown in Section 46, using a value close to the minimum number of
walkers, meaning twice the number of parameters, is generally a good choice.
Furthermore, we suggest to increase the number of walkers by a multiplicative
factor equal to the number of well separated modes (e.g. four times the number
of dimensions in a bimodal density). Other cases in which increasing the
number of walkers can improve the sampling efficiency include target
distributions with strong non-linear correlations between their parameters.

Regarding the initial positions of the walkers, we found that we can reduce
the length of the burn-in phase by initialising the walkers from a tight
sphere (i.e. Normal distribution with a very small variance) close to the
_Maximum a Posteriori_ (MAP) estimate. In high dimensional problems, the MAP
estimate will not reside in the typical set and the burn-in phase might be
longer. We found that the tight sphere initialisation is still an efficient
strategy compared to a more dispersed initialisation [50]. Other approaches
include initialising the walkers by sampling from the prior distribution or
the _Laplace approximation_ of the posterior distribution. In multimodal
cases, a prior initialisation is usually a better choice. A brief simulated
annealing phase can also be very efficient, particularly in cases with many
well separated modes.

Recent work on the No U-Turn Sampler [84] has attempted to reduce the hand-
tuning requirements of Hamiltonian Monte Carlo [15] using the dual averaging
scheme of [136]. In order to achieve a similar result, we employed the much
simpler stochastic approximation method of [148] to tune the initial length
scale $\mu$. The Affine Invariant Ensemble Sampler [72] and the Differential
Evolution MCMC [169] use an ensemble of walkers to perform Metropolis updates.
Our method differs by using the information from the ensemble to perform Slice
Sampling updates. So why does ESS perform better, as demonstrated, compared to
those other methods? The answer lies in the locally adaptive and non-rejection
nature of the algorithm (i.e. stepping out and shrinking) that enables both
efficient exploration of non-linear correlations and large steps in parameter
space (e.g. using the global move)888Indeed, large steps like the ones in the
50–dimensional Gaussian Mixture example would not have been possible without
the non-rejection aspect of the method as most attempts to jump to the other
mode would have missed it using Metropolis updates..

For all numerical benchmarks in this paper we used the publicly available,
open source Python implementation of Ensemble Slice Sampling called zeus999The
code is available at https://github.com/minaskar/zeus. [96].

#### 48 Conclusion

We have presented Ensemble Slice Sampling (ESS), an extension of Standard
Slice Sampling that eliminates the latter’s dependence on the initial value of
the length scale hyperparameter and augments its capacity to sample
efficiently and in parallel from highly correlated and strongly multimodal
distributions.

In this paper we have compared Ensemble Slice Sampling with the optimally-
tuned Metropolis and Standard Slice Sampling algorithms. We found that, due to
its affine invariance, Ensemble Slice Sampling generally converges faster to
the target distribution and generates chains of significantly lower
autocorrelation. In particular, we found that in the case of AR(1), Ensemble
Slice Sampling generates an order of magnitude more independent samples per
evaluation of the probability density than Metropolis and Standard Slice
Sampling. Similarly, in the case of the correlated funnel distribution,
Ensemble Slice Sampling outperforms Standard Slice Sampling by an order of
magnitude in terms of efficiency. Furthermore, in this case, Metropolis-based
proposals fail to converge at all, demonstrating that a single Metropolis
proposal scale is often not sufficient.

When compared to state-of-the-art ensemble methods (i.e. AIES, DEMC) Ensemble
Slice Sampling outperforms them by $1-2$ orders of magnitude in terms of
efficiency for target distributions with non-linear correlations (e.g. the
Ring and Gaussian shells distributions). In the real world example of
hierarchical Gaussian process regression, ESS’s efficiency is again superior
by $1-2$ orders of magnitude. Furthermore, in the Bayesian object detection
example ESS achieved higher accuracy compared to AIES and DEMC. Finally, in
the strongly multimodal case of the Gaussian Mixture, ESS outperformed all
other methods (i.e. SMC, AIES, DEMC, KM) and was the only sampler able to
produce reliable results in $50$ dimensions.

The consistent high efficiency of the algorithm across a broad range of
different problems along with its parallel, black-box and gradient-free
nature, renders Ensemble Slice Sampling ideal for use in scientific fields
such as physics, astrophysics and cosmology, which are dominated by a wide
range of computationally expensive and almost always non-differentiable
models. The method is flexible and can be extended further using for example
tempered transitions [89] or subspace sampling [182].

#### 49 Appendix: Estimating the Effective Sample Size

Assuming that the computational bottleneck of a MCMC analysis is the
evaluation of the probability density function, which is usually a valid
assumption in scientific applications, the _efficiency_ can be formally
defined as the ratio of the Effective Sample Size $N_{\rm Eff}$ to the total
number of probability evaluations for a given chain.

The $N_{\rm Eff}$ quantifies the number of effectively independent samples of
a chain, and it is defined as

$N_{\rm Eff}=\frac{n}{\text{IAT}}\,,$ (10.24)

where $n$ is the actual number of samples in the chain, and IAT is the
integrated autocorrelation time. The latter describes the number of steps that
the sampler needs to do in order to forget where it started and it is defined
as

$\text{IAT}=1+2\sum_{k=1}^{\infty}\rho(k)\,,$ (10.25)

where $\rho(k)$ is the normalised autocorrelation function at lag $k$. In
practice, we truncate the above summation in order to remove noise from the
estimate [164].

Given a chain $X(k)$ with $k=1,2,...,n$ the normalised autocorrelation
function $\hat{\rho}(k)$ at lag $k$ is estimated as

$\hat{\rho}(k)=\frac{\hat{c}(k)}{\hat{c}(0)}\,,$ (10.26)

where

$\hat{c}(k)=\frac{1}{n-k}\sum_{m=1}^{n-k}\big{[}X(k+m)-\bar{X}\big{]}\big{[}X(m)-\bar{X}\big{]}\,,$
(10.27)

and $\bar{X}$ is the mean of the samples.

In the case of ensemble methods, the IAT of an ensemble of chains is computed
by first concatenating the chain from each walker into a single long chain. We
found this estimator has lower variance than the [72] estimator and the [48]
estimator.

### 11 Zeus

This chapter presents zeus which is the main contribution introduced in the
paper titled zeus: A Python implementation of Ensemble Slice Sampling for
efficient Bayesian parameter inference that was published in the journal
Monthly Notices of the Royal Astronomical Society in 2021 [96]. The content of
the chapter is almost identical to that included in the aforementioned
publication with the exception of minor text and figure formatting
differences.

We introduce zeus, a well-tested Python implementation of the Ensemble Slice
Sampling (ESS) method for Bayesian parameter inference. ESS is a novel Markov
chain Monte Carlo (MCMC) algorithm specifically designed to tackle the
computational challenges posed by modern astronomical and cosmological
analyses. In particular, the method requires only minimal hand–tuning of $1-2$
hyper-parameters that are often trivial to set; its performance is insensitive
to linear correlations and it can scale up to 1000s of CPUs without any extra
effort. Furthermore, its locally adaptive nature allows to sample efficiently
even when strong non-linear correlations are present. Lastly, the method
achieves a high performance even in strongly multimodal distributions in high
dimensions. Compared to emcee, a popular MCMC sampler, zeus performs $9$ and
$29$ times better in a cosmological and an exoplanet application respectively.

#### 50 Introduction

Over the past few decades the volume of astronomical and cosmological data has
increased substantially. In response to that, a variety of astrophysical
models have been developed to explain the plethora of observations. Markov
chain Monte Carlo (MCMC) has been established as the standard procedure of
inferring the model parameters subject to the available data in a Bayesian
framework. Within the Bayesian context, the object that quantifies the
probability distribution of the model parameters $\theta$ given the data $D$
and model $\mathcal{M}$ is the posterior distribution
$\mathcal{P}(\theta)\equiv P(\theta|D,\mathcal{M})$ which is defined using
Bayes’s theorem:

$\mathcal{P}(\theta)=\frac{\mathcal{L}(\theta)\pi(\theta)}{\mathcal{Z}},$
(11.1)

where $\mathcal{L}(\theta)\equiv P(D|\theta,\mathcal{M})$ is the likelihood
function, $\pi(\theta)\equiv P(\theta|\mathcal{M})$ is the prior distribution
of the model parameters $\theta$, and $\mathcal{Z}\equiv P(D|\mathcal{M})$ is
the, so called, Bayesian model evidence or marginal likelihood and in this
context can be treated as a simple normalisation constant.

MCMC does not in general require knowing the value of the model evidence and
it only depends on the ability to evaluate the unnormalised posterior
distribution for arbitrary values of $\theta$. MCMC methods can then be used
to generate (Markov) chains of samples from the posterior distribution. Those
samples can be used to calculate integrals (e.g. parameter uncertainties,
marginal distributions etc.) that are paramount for modern astronomical and
cosmological analyses.

The most commonly used MCMC methods are variants of the Metropolis-Hastings
(MH) algorithm [122, 80]. MH consists of two steps. First, given the last
sample in the chain, a new sample is proposed and then the Metropolis
criterion determines whether or not that new sample should be accepted and
thus added to the chain. The resulting chain is Markovian in the sense that
each sample is proposed based only on the previous sample. The purpose of the
Metropolis acceptance criterion is to bias the chain so that the time spent in
a region of the parameter space would be proportional to the posterior
probability in that region. In other words, the stationary distribution of the
Markov chain is the target distribution i.e. the posterior distribution. For a
detailed introduction to MCMC methods we direct the reader to [117] and for an
intuitive introduction to Bayesian inference to [92].

Arguably, the most difficult part of the MH algorithm is the proposal step.
There are many ways of choosing a new sample and the efficiency of the method
depends on this choice. By far the simplest one is the use of a normal
(Gaussian) distribution, centred around the previous sample to generate the
new proposed sample. The resulting method is often called Random Walk
Metropolis algorithm and its performance is highly sensitive to the $n(n+1)/2$
elements that form its covariance matrix. Those elements generally need to be
chosen a priori or be adaptively tuned. More efficient methods utilise the
gradient of the target distribution [14] or an ensemble of parallel and
communicating chains [69, 169, 170, 72].

Out of the methods mentioned in the previous paragraph we will focus our
attention on the last one, the ensemble or population MCMC variety. The reason
is simple: the Random Walk Metropolis algorithm requires a great amount of
tuning (or a priori knowledge) for it to perform efficiently and even then
there is no guarantee that the proposal covariance matrix is optimal for the
whole parameter space. On the other hand, gradient based methods, although
very powerful, are in general unsuitable for astronomical applications in
which the models that are used are almost always not differentiable.

One benefit of ensemble MCMC over its alternatives is that the ensemble of
parallel chains (also known as walkers) collectively sample the posterior,
thus information about their distribution can be shared and used to make
better educated proposals. Other advantages include the lack of hand-tuning of
hyper-parameters and their capacity for parallel implementation. For the
aforementioned reasons, ensemble MCMC methods have dominated astronomical
analyses. The most common ones are affine–invariant ensemble sampling (AIES)
[72] and differential evolution MCMC (DEMC) [169, 170], both implemented in
the popular Python package emcee [50, 49].

In this paper we introduce zeus, a stable and well-tested Python
implementation of Ensemble Slice Sampling (ESS) [95]. ESS is a method based on
the ensemble MCMC paradigm, with the crucial difference being that its
proposals are performed via Slice Sampling updates [133] instead of
Metropolis-Hastings ones. As we will thoroughly demonstrate in Section 52,
this subtle difference leads to substantial improvements in terms of sampling
efficiency and robustness. zeus is a user-friendly tool that does not require
any hand-tuning or preliminary runs and can scale up to 1000s of CPUs without
any extra effort from the user.

zeus has been used in various astronomical and cosmological analyses,
including cosmological tests of gravity [168], relativistic effects and
primordial non-Gaussianity [183], 21cm intensity mapping [176], and has been
implemented as part of the CosmoSIS package [188].

zeus is open source software that is publicly available at
https://github.com/minaskar/zeus under the GPL-3 Licence. Detailed
documentation and examples on how to get started are available at
https://zeus-mcmc.readthedocs.io.

#### 51 Ensemble Slice Sampling

zeus is a Python implementation of the Ensemble Slice Sampling (ESS) method
presented in [95]. Here we will provide a high-level description of the method
and will refer to the accompanying paper for more details about the underlying
algorithmic structure and mathematics.

ESS combines the ensemble MCMC paradigm with slice sampling. Since the use of
slice sampling in astronomical parameter inference is rare we will start by
explaining its function and how it differs from Metropolis updates. Then we
will move on to discuss how it can be efficiently combined with ensemble MCMC.

##### 51.1 Slice sampling

Slice sampling is based on the idea that sampling from a distribution with
density $P(x)$ is equivalent to uniform sampling from the area under the plot
of $f(x)\propto P(x)$. To this end, we introduce an auxiliary variable $y$,
called height, such that the joint distribution $P(x,y)$ is uniform over the
region $U=\\{(x,y):0<y<f(x)\\}$. To sample from the marginal distribution
$P(x)$, we first sample from $P(x,y)$ and then we marginalise by dropping the
$y$ value of each sample.

In order to generate samples from $P(x,y)$ we utilise the following scheme
[133]:

  1. 1.

Given the current state $x_{0}$, draw $y_{0}$ uniformly from $(0,f(x_{0}))$.

  2. 2.

Find an interval $I=(L,R)$ that contains all, or at least part, of the slice
$s=\\{x:y_{0}<f(x)\\}$.

  3. 3.

Draw the new sample $x_{1}$ uniformly from $I\cap S$.

Figure 11.1: Illustration of the univariate slice sampling update. Given the
current sample $x_{0}$, a value $y_{0}$ is uniformly sampled along the
vertical slice $(0,f(x_{0}))$ (dashed line) thus defining the initial point
(blue). An interval $(L,R)$ is uniformly positioned horizontally around
$(x_{0},y_{0})$ and it is expanded in steps of size $R-L$ until both its ends
are outside the slice. The new sample is generated by repeatedly sampling
(uniformly) from the interval $(L^{\prime},R^{\prime})$ until a sample (green
star) is found inside the slice. Samples outside of the slice (red star) are
rejected and they are instead used to shrink $(L^{\prime},R^{\prime})$.

To construct the interval $I$ (step ii), [133] introduced the stepping-out
procedure that works by randomly positioning an interval of length $\mu$
around the sample $x_{0}$ (i.e. blue dot in Figure 11.1) and then expanding it
in steps of size $\mu$ until both its ends (i.e. $L^{\prime}$ and
$R^{\prime}$) are outside the slice. To obtain $x_{1}$ (i.e. green star in
Figure 11.1) we then use the shrinking procedure in which candidates are
sampled uniformly from $I$ until a point inside the slice $S$ is found.
Samples outside of the slice are used to shrink the interval $I$. The two
procedures are shown in Figure 11.1.

The length scale $\mu$ is the only free hyperparameter of slice sampling and
although its choice can reduce or increase the computational cost of the
method it generally does not affect its mixing properties (e.g. convergence
rate, autocorrelation time, etc.). zeus utilises a stochastic optimization
algorithm similar to [171] and based on the [148] optimisation scheme in order
to tune $\mu$ to its optimal value (see Section 3.1 of [95] for more details).

It is important to note here that for multimodal target distributions there is
no guarantee that the approximate slice would cross any of the other modes. In
particular, if the initial estimate of the length scale $\mu$ is low then the
probability of missing the other peaks, assuming that they are located far
away, is also low. As we will show in Section 52, unlike simple slice
sampling, ESS and thus zeus does not suffer from this effect.

##### 51.2 Walkers, moves and parallelism

Figure 11.2: The figure illustrates the differential move in the context of
Ensemble Slice Sampling. The walker $X_{k}$ to be updated is shown in red. Two
walkers, $X_{l}$ and $X_{m}$, (blue) are uniformly selected from the
complementary ensemble (grey). The approximate slice (dotted line) is
constructed parallel to the two walkers $X_{l}$ and $X_{m}$ using the
stepping-out procedure. The new position $Y$ (green) of $X_{k}$ is sampled
using the shrinking procedure along the approximate slice.

The slice sampling update described in the previous paragraphs is a univariate
update scheme. For it to be used to sample from multivariate target
distributions it needs to be generalised accordingly. Perhaps the simplest
such generalisation in a multivariate setting is the use of slice sampling to
sample along each coordinate axis in turn (i.e. component-wise slice sampling)
or to sample along randomly selected directions in parameter space [117].
Although valid, both of these approaches are unsuitable in cases of correlated
parameters in which the proper choice of direction can substantially
accelerate mixing.

To address this issue, [171] proposed to orthogonalise the parameter space
using the sample covariance, thus getting rid of linear correlations between
parameters. We will instead follow a different, perhaps more flexible,
approach to construct an efficient slice sampler. Our aim is to utilise an
ensemble of parallel chains/walkers that can exchange information about the
covariance structure of the target distribution and thus by-pass the
difficulties posed by correlations.

As hinted in the introduction, the ensemble of walkers collectively sample the
target distribution and thus their positions encode information about the
correlations between the parameters. One way to take advantage of this
information is to use it to construct direction vectors along which slice
sampling can take place. Many moves that generate direction vectors from the
complementary ensemble are possible. zeus offers a collection of them,
including some that utilise clustering algorithms and density estimation
methods. As we will show in Section 52, such moves can help accelerate
sampling in difficult cases such as strongly multimodal distributions. Any
distribution of the complementary ensemble can be used as a valid proposal to
generate such direction vectors and zeus offers a highly flexible interface
for the user to define such a move or choose one (or a mixture) from the ones
that are already implemented and tested. Here is a list of the currently
implemented moves in zeus:

  * $\bullet$

Differential move: This is the default move used by zeus and shown in Figure
11.2. Using the differential move, Ensemble Slice Sampling updates the
position of each walker in the ensemble by slice sampling along a direction
defined by the difference between two uniformly selected walkers from the rest
of the ensemble (i.e. the complementary ensemble).

  * $\bullet$

Gaussian move: The Gaussian move samples the direction vectors along which
slice sampling is performed from a normal distribution that shares the same
covariance structure as the complementary ensemble. This approach is very
efficient in cases in which the target distribution is close to normal.

  * $\bullet$

Global move: The Global move utilises a Dirichlet Process Gaussian Mixture to
fit the complementary ensemble and proposes directions along different peaks
of the target distribution in cases of strong multi-modality.

  * $\bullet$

KDE move: The KDE move samples the direction vectors from a Gaussian Kernel
Density Estimate of the complementary ensemble. This can be useful in cases of
highly non-Gaussian target distributions.

  * $\bullet$

Random move: The Random move performs slice sampling along isotropic
directions. This is equivalent of standard multivariate slice sampling and it
is mostly offered for testing purposes as it cannot handle correlations
efficiently.

For more information on how those moves work as well as a comparison of the
Differential, Gaussian and Global moves we direct the interested reader to