# The stepping-stone sampling algorithm for calculating the evidence of
gravitational wave models

Patricio Maturana Russel1, Renate Meyer1, John Veitch2 and Nelson
Christensen3,4 1 Department of Statistics, University of Auckland, Auckland
1142, New Zealand
2 Institute for Gravitational Research, School of Physics and Astronomy,
University of Glasgow, Glasgow G12 8QQ, United Kingdom
3 ARTEMIS, Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, CS
34229, F-06304 Nice Cedex 4, France
4 Physics and Astronomy, Carleton College, Northfield, MN 55057, USA

###### Abstract

Bayesian statistical inference has become increasingly important for the
analysis of observations from the Advanced LIGO and Advanced Virgo
gravitational-wave detectors. To this end, iterative simulation techniques, in
particular nested sampling and parallel tempering, have been implemented in
the software library LALInference to sample from the posterior distribution of
waveform parameters of compact binary coalescence events. Nested sampling was
mainly developed to calculate the marginal likelihood of a model but can
produce posterior samples as a by-product. Thermodynamic integration is
employed to calculate the evidence using samples generated by parallel
tempering but has been found to be computationally demanding. Here we propose
the stepping-stone sampling algorithm, originally proposed by Xie et al.
(2011) in phylogenetics and a special case of path sampling, as an alternative
to thermodynamic integration. The stepping-stone sampling algorithm is also
based on samples from the power posteriors of parallel tempering but has
superior performance as fewer temperature steps and thus computational
resources are needed to achieve the same accuracy. We demonstrate its
performance and computational costs in comparison to thermodynamic integration
and nested sampling in a simulation study and a case study of computing the
marginal likelihood of a binary black hole signal model applied to simulated
data from the Advanced LIGO and Advanced Virgo gravitational wave detectors.
To deal with the inadequate methods currently employed to estimate the
standard errors of evidence estimates based on power posterior techniques, we
propose a novel block bootstrap approach and show its potential in our
simulation study and LIGO application.

thermodynamic integration, stepping-stone sampling, nested sampling, model
selection, gravitational waves

###### pacs:

04.30.-w, 02.50.-r, 05.45.Tp, 97.60.Bw

##  I Introduction

It has now been two decades since Bayesian parameter estimation routines were
first introduced for studies in astrophysics Saha and Williams (1994),
gravitational waves Christensen and Meyer (1998), and cosmology Christensen
_et al._ (2001); Knox  _et al._ (2001). Bayesian parameter estimation routines
have become extremely important for these disciplines, and their use is
ubiquitous Sharma (2017). Recent dramatic observations with, for example, the
cosmic microwave background and gravitational waves, have been used with
Bayesian parameter estimation methods to significantly push our knowledge of
the universe and its history Hinshaw  _et al._ (2013); Planck Collaboration
_et al._ (2014, 2016); Abbott  _et al._ (2016a, b); Meyer and Christensen ;
Abbott  _et al._ (2017a). Advances in computer power, coupled with new and
innovative Bayesian parameter estimation techniques, continue to push the
applicability and importance of Bayesian methods Zhu  _et al._ (2017).

The importance of accurate parameter estimation calculations was dramatically
displayed with the observations of gravitational waves and gamma rays from the
binary neutron star merger GW170817 and GRB 170817A Abbott  _et al._ (2017a);
Goldstein  _et al._ (2017). Using the data from the two Advanced LIGO
detectors Aasi  _et al._ (2015) and the Advanced Virgo detector Acernese  _et
al._ (2015) an initial sky-map and distance estimates from the gravitational
wave data was released five hours after the merger Abbott  _et al._ (2017b)
using a specially designed method for sky position estimation Singer and Price
(2016). A little over 11 hours after the gravitational wave - gamma ray event
refined estimates were released based on the first comprehensive parameter
estimation Veitch  _et al._ (2015), giving a more accurate estimated of the
sky position and distance to the source Abbott  _et al._ (2017b). The
parameter estimation calculations allowed astronomers to identify the source,
providing for electromagnetic observations that yielded a plethora of
astrophysical information, including the observation of a kilonova Abbott  _et
al._ (2017b). Parameter estimation of gravitational wave models will continue
to be significant for multimessenger astronomy.

Also of critical importance is the ability to conduct model comparison and
parameter estimation studies with the gravitational wave signals. For example,
Bayesian parameter estimation methods were used to decipher the physical
characteristics of the observed gravitational wave events, such as the first
observed binary black hole merger GW150914 Abbott  _et al._ (2016a, b, c, d),
and the binary neutron star merger GW170817 Abbott  _et al._ (2017a, 2018a).
Similarly, model comparisons have been conducted in a number of ways using the
data from the detected gravitational waves signals. This includes tests of
general relativity Abbott  _et al._ (2016e), neutron star equation of state
studies Abbott  _et al._ (2017a, 2018b), constraining tidal instabilities in
binary neutron star mergers Weinberg (2018), and the search for a stochastic
gravitational wave background from binary black hole mergers over the history
of the universe Smith and Thrane (2018). When Advanced LIGO and Advanced Virgo
made the first observation of a binary black merger using the data from three
detectors it provided an opportunity to conduct a model comparison test as to
whether the polarization of the gravitational waves was consistent with
general relativity, or other theories of gravity; general relativity succeeded
in this important model comparison Abbott  _et al._ (2017c). Advanced LIGO
data was also used to search for a stochastic gravitational wave background as
described by general relativity or alternative theories of gravity, and model
comparison was integral to this study Callister  _et al._ (2017); Abbott  _et
al._ (2018c). Methods that improve the calculation of the evidence, the
marginal likelihood of a model, would be well-received in the gravitational
wave community, and would certainly be of use in other areas of astrophysics
and cosmology Sharma (2017). Equally important is an accurate estimation of
the associated standard error. Here we introduce the moving block bootstrap
(MBB) that accounts for the autocorrelation between the samples, and provides
a more accurate estimate than the standard bootstrap method.

Presented in this paper is the stepping-stone sampling (SS) algorithm that
provides an improvement of the evidence estimator for Bayesian model selection
and the MBB for computing its standard error. The practicality of the SS
algorithm for calculating the evidence of gravitational wave models will be
demonstrated. The SS algorithm could be a further advancement for model
selection for gravitational waves data analysis, as well as for other
applications in astrophysics and astronomy.

The paper is structured as follows. In Section II we review nested sampling,
thermodynamic integration and introduce the stepping-stone algorithm for
computing the evidence of Bayesian model selection. In Section III we
introduce the moving block bootstrap for calculating the Monte Carlo standard
error of the evidence estimates. In Section IV we show the enhanced
performance of the SS algorithm over thermodynamic integration in a simulation
study. The different algorithms are then applied to simulated LIGO-Virgo
gravitational wave data in Section V. Their results are contrasted and the
benefits of the SS algorithm and the moving black bootstrap for standard error
estimation become evident. A summary discussion is given in Section VI.

##  II Computation of Marginal Likelihood

The evidence or marginal likelihood of a model $M$ is a multi-dimensional
integral defined as

$\displaystyle
z=\int_{\Theta}L(\bm{X}|\bm{\theta},M)\pi(\bm{\theta}|M)\text{d}\bm{\theta},$
(1)

where $\bm{\theta}\in\Theta$ denotes the parameter vector, $\bm{X}$ the
dataset, $L(\bm{X}|\bm{\theta},M)$ the likelihood function, and
$\pi(\bm{\theta}|M)$ the prior density, assumed to be proper, i.e.
$\int_{\Theta}\pi(\bm{\theta}|M)\text{d}\bm{\theta}=1$.

In general, this integral (1) has no analytical solution and must be estimated
using numerical methods. Importance sampling techniques, in particular the
arithmetic mean (AM) and harmonic mean (HM) methods, provide the simplest way
of estimating it (Newton and Raftery, 1994). Let
$\bm{\theta}_{i},i=1,\ldots,n$ be samples from the prior, the AM estimator is
an average of corresponding $n$ likelihood values:

$\displaystyle\widehat{z}_{AM}=\frac{1}{n}\sum_{i=1}^{n}L(\bm{X}|\bm{\theta}_{i},M).$
(2)

In general, high-likelihood areas are very small and constitute a small
fraction of the prior. Therefore, unless $n$ is very large, the sample will
not adequately represent these areas and yield a poor estimate. The HM
estimator is based on samples $\bm{\theta}_{i},i=1,\ldots,n$ drawn from the
posterior:

$\displaystyle\widehat{z}_{HM}=\left(\frac{1}{n}\sum_{i=1}^{n}\frac{1}{L(\bm{X}|\bm{\theta}_{i},M)}\right)^{-1}~{}.$
(3)

This is the harmonic mean of likelihood values.

The AM and HM estimators are not recommended because they produce unreliable
estimates of the evidence, even though they are easily calculated. In this
context, more complex approaches have been proposed, such as power posterior
methods Xie  _et al._ (2011); Lartillot and Philippe (2006); Friel and Pettitt
(2008); Neal (2001). These methods rely on a set of transitional distributions
which connect the prior and the posterior, reminiscent of simulated annealing.
The geometric path is the most popular scheme used to connect these
distributions and defines the power posterior density as

$\displaystyle
p_{\beta}(\bm{\theta}|\bm{X},M)=\frac{L(\bm{X}|\bm{\theta},M)^{\beta}\pi(\bm{\theta}|M)}{z_{\beta}},$
(4)

for the inverse temperature $0\leq\beta\leq 1$, where $z_{\beta}$ is the
normalizing constant, which is defined as
$\int_{\Theta}L(\bm{X}|\bm{\theta},M)^{\beta}\pi(\bm{\theta}|M)\text{d}\bm{\theta}$.
Note that the power posterior density turns into the prior and posterior for
$\beta=0$ and $\beta=1$, respectively.

Methods that make use of samples from the power posteriors are much more
accurate than HM as has been widely documented (Lartillot and Philippe, 2006;
Newton and Raftery, 1994; Xie  _et al._ , 2011), particularly in high
dimensional problems. Among these methods, thermodynamic integration (TI)
(Lartillot and Philippe, 2006) is a popular method to estimate the evidence of
gravitational wave (GW) models, showing in general good performance. Another
method, widely applied in other fields such as phylogenetics is the SS
algorithm (Xie  _et al._ , 2011). As this method can provide many advantages
over the TI estimate, it is important to explore the performance of the SS
estimator for GW models as to the best of our knowledge, the SS algorithm has
not been used for evidence calculation in this context.

One of the drawbacks of power posterior methods is the significant
computational cost required to produce a single evidence estimate as multiple
Markov chains have to be run, one for each temperature. Fortunately, since
parallel tempering is commonly used in GW parameter estimation, the samples at
different temperatures are available and can be recycled in order to use these
methods.

However, as has been noticed in Veitch  _et al._ (2015), TI might require a
larger number of temperatures than the one needed for parameter estimation in
order to achieve accurate estimates. Note that the samples of chains at
temperatures $T>1$ ($\beta<1$) are only used to aid the mixing of the chain at
$T=\beta=1$ whose stationary distribution is the posterior, and are therefore
discarded from the inference process. In this context, the SS algorithm seems
very promising since it requires fewer temperature steps than TI to provide
accurate evidence estimates as we will show in section IV.

Another method to estimate the evidence, not based on power posteriors, is
nested sampling (NS) (Skilling, 2006; Veitch and Vecchio, 2010). This Bayesian
algorithm has been successfully applied in diverse fields, such as astronomy
(Brewer and Donovan, 2015), cosmology (Feroz  _et al._ , 2009), engineering
(Henderson  _et al._ , 2017) and phylogenetics (Maturana Russel  _et al._ ,
2018; Maturana Russel, 2018). To estimate the evidence of GW models, NS has
been implemented in the software package LALInference (Veitch  _et al._ ,
2015). The method has the unique property of yielding an estimation of the
uncertainty associated with the evidence estimate in a single run (however,
only for independent samples).

Alternatively, instead of estimating the evidence for each model being tested,
a trans-dimensional Reversible Jump Markov chain Monte Carlo (RJMCMC; Green,
1995; Umstätter  _et al._ , 2005) method can be used in order to explore the
joint space of all models. Then the probability for each model can be
calculated simply by calculating the relative frequency of visits to each
model by the Markov chain. However, this exploration depends on tuning
parameters which can be difficult to specify, leading to poor mixing of the
Markov chain and subsequently to large statistical errors associated with the
evidence estimates (Cornish and Littenberg, 2015).

Below we describe TI, SS and NS in more detail before comparing their
performance in sections IV and V.

###  II.1 Thermodynamic Integration

Thermodynamic integration or the more general path sampling (Gelman and Meng,
1998) make use of an auxiliary variable $\beta$, $0\leq\beta\leq 1$, to define
transitional distributions, namely the power posterior distributions defined
in (4) in the case of TI, that provide a path from the prior ($\beta=0$) to
the posterior distribution ($\beta=1$). By explicitly denoting the evidence
$z_{\beta}$ as a function of $\beta$ by

$\displaystyle
z(\bm{X}|\beta)=\int_{\Theta}L(\bm{X}|\bm{\theta},M)^{\beta}\pi(\bm{\theta}|M)\text{d}\bm{\theta},$
(5)

the log marginal likelihood has the representation as the integral over the
1-dimensional parameter $\beta$ of half the mean deviance where the
expectation is taken with respect to the power posterior:

$\log(z)=\log\left(\frac{z(\bm{X}|\beta=1)}{z(\bm{X}|\beta=0)}\right)=\int_{0}^{1}E_{\beta}\left[\log(p(\bm{X}|\bm{\theta},M)\right]\text{d}\beta.$
(6)

Representation (6) follows by integration from

$\displaystyle\frac{\partial}{\partial\beta}\log(z(\bm{X}|\beta))=\frac{1}{z(\bm{X}|\beta)}\frac{\partial}{\partial\beta}z(\bm{X}|\beta)$
$\displaystyle=\frac{1}{z(\bm{X}|\beta)}\frac{\partial}{\partial\beta}\int_{\Theta}L(\bm{X}|\bm{\theta},M)^{\beta}\pi(\bm{\theta}|M)\text{d}\bm{\theta}$
$\displaystyle=\frac{1}{z(\bm{X}|\beta)}\int_{\Theta}L(\bm{X}|\bm{\theta},M)^{\beta}\log(L(\bm{X}|\bm{\theta},M))\pi(\bm{\theta}|M)\text{d}\bm{\theta}$
$\displaystyle=\int_{\Theta}\frac{L(\bm{X}|\bm{\theta},M)^{\beta}\pi(\bm{\theta}|M)}{z_{\beta}}\log(L(\bm{X}|\bm{\theta},M))\text{d}\bm{\theta}$
$\displaystyle=E_{\beta}\left[\log(L(\bm{X}|\bm{\theta},M)\right].$

The samples from the parallel tempered chains for different values of $\beta$
provide samples from the power posteriors and the expectation
$E_{\beta}\left[\log(L(\bm{X}|\bm{\theta},M)\right]$ is then estimated by the
sample average. The integral in equation (6) is then approximated by numerical
integration, e.g. using the trapezoidal or Simpson’s rule.

###  II.2 Stepping-stone Sampling Algorithm

Stepping-stone sampling is another method to estimate the marginal likelihood.
It has been widely used by the phylogenetic community where it was proposed by
Xie  _et al._ (2011). SS works basically by mixing elements from importance
sampling and simulated annealing methods. This method relies on the same
sampling scheme required by TI. Therefore, its implementation in any software
package where TI or parallel tempering has already been implemented should be
straightforward. SS has the advantage of requiring fewer path steps than TI to
accurately estimate the marginal likelihood and yielding a less-biased
estimator as demonstrated in section IV.

The marginal likelihood can be seen as the ratio $z=z_{1}/z_{0}$, where
$z_{0}=1$ since the prior is assumed to be proper. The direct calculation of
this ratio via importance sampling is not reliable because the distributions
involved in the numerator and denominator (posterior and prior, respectively)
are, in general, quite different. To solve this problem, SS expands this ratio
in a telescope product of $K$ ratios of normalizing constants of the
transitional distributions Neal (1993), that is

$\displaystyle
z=\frac{z_{1}}{z_{0}}=\frac{z_{\beta_{1}}}{z_{\beta_{0}}}\frac{z_{\beta_{2}}}{z_{\beta_{1}}}\dots\frac{z_{\beta_{K-2}}}{z_{\beta_{K-3}}}\frac{z_{\beta_{K-1}}}{z_{\beta_{K-2}}}=\prod_{k=1}^{K-1}\frac{z_{\beta_{k}}}{z_{\beta_{k-1}}}=\prod_{k=1}^{K-1}r_{k},$

for $\beta_{0}=0<\beta_{1}<\dots<\beta_{K-2}<\beta_{K-1}=~{}1$, being the
sequence of inverse temperatures, where $r_{k}=z_{\beta_{k}}/z_{\beta_{k-1}}$.
These individual intermittent ratios can be estimated with higher accuracy
than $\frac{z_{1}}{z_{0}}$ because the distributions in the numerator and
denominator are generally quite similar when using a reasonable number of
temperatures $K$. In this situation the importance sampling method works well.

SS estimates each ratio $r_{k}$ by importance sampling using $p_{\beta_{k-1}}$
as importance sampling distribution. This is a suitable distribution because
it has heavier tails than $p_{\beta_{k}}$ which leads to an efficient estimate
of $r_{k}$. In this manner, it avoids estimating from the posterior
distribution, making it slightly less expensive computationally than TI for
the same number of path steps. The estimation of each ratio is based on the
identity

$\displaystyle r_{k}=\frac{z_{\beta_{k}}}{z_{\beta_{k-1}}}$
$\displaystyle=\int_{\Theta}\frac{L(\bm{X}|\bm{\theta},M)^{\beta_{k\quad}}}{L(\bm{X}|\bm{\theta},M)^{\beta_{k-1}}}\>p_{\beta_{k-1}}(\bm{\theta}|\bm{X},M)\text{d}\bm{\theta},$

which is estimated by its unbiased Monte Carlo estimator

$\displaystyle\widehat{r}_{k}=\frac{1}{n}\sum_{i=1}^{n}L(\bm{X}|\bm{\theta}_{\\!\beta_{k-1}}^{i},M)^{\beta_{k}-\beta_{k-1}},$

where
$\bm{\theta}_{\\!\beta_{k-1}}^{1},\dots,\bm{\theta}_{\\!\beta_{k-1}}^{n}$ are
drawn from $p_{\beta_{k-1}}$ with $k=1,\dots,K-1$.

Therefore, the SS estimate of the marginal likelihood is defined as

$\displaystyle\widehat{z}$
$\displaystyle=\prod_{k=1}^{K-1}\frac{1}{n}\sum_{i=1}^{n}L(\bm{X}|\bm{\theta}_{\\!\beta_{k-1}}^{i},M)^{\beta_{k}-\beta_{k-1}},$

with log-version

$\displaystyle\log\widehat{z}$
$\displaystyle=\sum_{k=1}^{K-1}\log\sum_{i=1}^{n}L(\bm{X}|\bm{\theta}_{\\!\beta_{k-1}}^{i},M)^{\beta_{k}-\beta_{k-1}}-(K-1)\log
n.$

Although $\widehat{z}$ is unbiased, the log transformation introduces a bias
which can be alleviated by increasing $K$ (Xie  _et al._ , 2011).

The performance of this method depends naturally on its specifications such as
the number of transitional distributions and number of samples from each of
them ($K$ and $n$, respectively). The dispersal of the $\beta$ values has also
a strong influence, even more so in TI (see Xie  _et al._ (2011) and our
simulation study below). Along these lines, Xie  _et al._ (2011) proposed to
spread the $\beta$ values according to the evenly spaced quantiles of a
Beta(0.3, 1) distribution. This distribution is right skewed, thereby putting
half of the $\beta$ values below 0.1 where most of the variability is found.

SS is closely related to annealed importance sampling (Neal, 2001). The latter
utilizes the same product of ratios, but instead of estimating each ratio
separately, it estimates the entire product via importance sampling, that is
the whole telescope product is evaluated multiple times and then these values
are averaged (Maturana Russel, 2017). For the particular case of $K=2$, that
is considering only the prior, both methods reduce to the arithmetic mean, and
for $n=1$, they are equivalent.

###  II.3 Nested Sampling

NS transforms the multidimensional integral defined in (1), by making use of a
property of positive random variables (see Maturana Russel  _et al._ (2018)
for more details), into a one-dimensional one that utilizes a function that
relates the prior with the likelihood as

$\displaystyle z=\int_{0}^{1}L(\xi)\text{d}\xi,$

where $L$ is the likelihood as a function of the prior volume $\xi$. This
function can be read as the proportion of prior volume $\xi$ with likelihood
values greater than $L(\xi)$.

This likelihood is a non-increasing function over the unit range. For a given
decreasing sequence of $\xi$-values and an increasing sequence of $L$-values,
the marginal likelihood can be estimated using, for instance, the trapezium
rule

$\displaystyle\widehat{z}_{NS}=\sum_{i=1}^{K}\dfrac{1}{2}(\xi_{i-1}-\xi_{i+1})L_{i},$

where $0<\xi_{K+1}<\xi_{K}<\cdots<\xi_{1}<\xi_{0}=1$.

NS explores the parameter space from the prior toward those areas of high
likelihood values over time. For this, a set of $N$ points, called live
points, are drawn independently from the prior. The point $\bm{\theta}_{1}$
with the lowest likelihood associated to these points is detected and the
latter is registered as $L_{1}$. Then, this point $\bm{\theta}_{1}$ is
replaced by a new one $\bm{\theta}^{*}$ drawn from the prior but restricted to
have a greater likelihood, that is $L(\bm{\theta}^{*})>L(\bm{\theta}_{1})$.
This procedure is repeated until a given stopping criterion is satisfied.
Thus, an increasing sequence of likelihood values $L_{1},\dots,L_{K}$ is
generated.

Even though the $\xi$-values cannot be measured precisely, the nature of this
algorithm allows them to be estimated. The $\xi$-sequence can be defined as

$\displaystyle\xi_{1}=u_{1},\>\xi_{2}=u_{2}\xi_{1},\dots,\>\xi_{K}=u_{K}\xi_{K-1},$

where $u_{i}\sim\text{Beta}(N,1)$. The geometric mean is the most common
method to estimate the $u$-values, which yields

$\displaystyle\xi_{i}=e^{-i/N}.$

The nature of NS algorithm also allows to estimate the standard error of the
$\log z$ estimate in a single run as

$\widehat{s.e.}_{\text{NS}}(\log z)=\sqrt{\dfrac{H}{N}},$ (7)

where $H$ is the negative entropy. However, this NS standard error estimate is
only valid if the samples are drawn independently. In practice though, the
samples will often be serially dependent because Metropolis-Hastings
algorithms are used for their generation. As an alternative to (7), for a
fixed sequence of likelihood values and multiple sequence of $\xi$-values,
generated from different $u\\!\sim\\!\text{Beta}(N,1)$ values, a distribution
of marginal likelihood estimates can be generated and subsequently the
uncertainty can be estimated.

##  III Estimation of the Monte Carlo Standard Error of the Evidence

The point estimate of the evidence is subject to random errors and therefore
we need to have a measure of the Monte Carlo standard error of the evidence
estimates. This is also important if we want to compare the performance of
different types of evidence estimates. In the NS case, the algorithm provides
direct ways of calculating its standard error from a single run as given in
(7). However, power posterior methods lack a reliable direct way of
calculating the standard error of the evidence. In Lartillot and Philippe
(2006) and Xie  _et al._ (2011), the authors proposed estimates which rely on
the independence of the samples in the Markov chains at different
temperatures, an assumption that is not met in general. Practitioners opt for
the standard procedure of repeating the analysis multiple times and then
calculating the standard error. This brute force technique can be very costly
and is in some cases computationally not viable. Alternatively, some estimate
the error internally in a single run, that is by re-sampling independently the
Markov chains in order to generate multiple evidence estimates. However, this
approach does not consider the potential autocorrelation in the samples,
leading to wrong estimates. Here, we propose the use of a block bootstrap
method for multivariate time series, which accounts for the autocorrelation
between the samples within a Markov chain at a fixed temperature and the
cross-correlation between parallel chains at different temperatures.

Bootstrap is a resampling procedure proposed by Efron (1979), initially for
independent variables and later generalized by several authors. An extension
for the case of time series was proposed in Kunsch (1989), which differs from
the original algorithm by allowing the sampling in blocks. The method is known
as moving block bootstrap, in short MBB. This allows to take into account the
presence of dependence in the data.

Let $X_{1},\dots,X_{n}$ be the observed values from a sequence of stationary
random variables, in our case, a Markov chain. Define the overlapping blocks
$B_{i}=(X_{i},,\dots,X_{i+\ell-1})$ of length $\ell$, for $1\leq i\leq
n-\ell+1$ and $1\leq\ell\leq n$, that is

$\displaystyle B_{1}$ $\displaystyle=(X_{1},X_{2},X_{3},\dots,X_{\ell})$
$\displaystyle B_{2}$ $\displaystyle=(X_{2},X_{3},X_{4},\dots,X_{\ell+1})$
$\displaystyle\>\>\vdots$ $\displaystyle B_{m}$
$\displaystyle=(X_{n-\ell+1},\dots,X_{n}),$

where $m=n-\ell+1$. MBB works by resampling randomly $b$ blocks (for didactic
reasons, suppose that $b=n/\ell$) and concatenating them in order to form a
set of bootstrap observations $X^{*}_{1},\dots,X^{*}_{n}$. For $\ell=1$, the
original bootstrap method for i.i.d. data is recovered. This procedure is
repeated as usual, generating the distribution of the statistic of interest,
in our case the marginal likelihood. In the general case that $n$ is not a
multiple of $\ell$, we can concatenate the random sample of $b$ block
bootstraps, where $b$ is $n/\ell$ rounded up, and discard the leftover points
$X^{*}_{n+1},\dots,X^{*}_{b\ell}$, such that the bootstrap observation set has
length $n$, as the original dataset.

Variants of this method can be found in Lahiri (2003), such as stationary
bootstrap, where the block length follows a geometric distribution;
nonoverlapping block bootstrap, which as its name say, considers
nonoverlapping blocks; and circular block bootstrap, which increases the
original dataset with the first $\ell-1$ observations in order to give equal
weights to all of them.

In the context of parallel tempering, in which case there are multiple Markov
chains, we need to generate the bootstrap observations using the same scheme
for all the chains. For instance, assuming equal chain lengths, a bootstrap
observation set for a Markov chain consisting in $(B_{6},B_{4},B_{2})$, is
replicated across the other chains. This procedure takes into account the
potential autocorrelation within the chains and the cross-correlation between
the chains due to the swaps in parallel tempering sampling. This is the
approach applied in our examples.

##  IV Simulation Study

We consider a simple Gaussian model used by Lartillot and Philippe (2006) to
test TI and compare it to the harmonic mean method. Here, it is used to
compare SS to TI. We also assess the error estimate via the MBB method and
compare it to the empirical calculation of the error. In addition, we study NS
performance for different sampling specifications.

The model is parametrized by a vector $\bm{x}=(x_{1},x_{2},\dots,x_{d})$ of
dimension $d$. The prior on $\bm{x}$ is a product of independent standard
normal distributions on each $x_{i}$, for $i=1,\dots,d$. The likelihood is

$\displaystyle L(\bm{x})=\prod_{i=1}^{d}e^{-\frac{x_{i}^{2}}{2v}},$

where $v$ is a parameter. Doing some calculations, it is easy to see that the
posterior distribution is given by a product of independent
$\text{N}(0,v/(1+v))$ distributions, and therefore, its marginal likelihood
has an analytical solution, which is $z=(v/(1+v))^{d/2}$. The power posterior
or transitional distributions are given by a product of independent
$\text{N}(0,v/(v+\beta))$ distributions. All the involved distributions are
Gaussians, so the sampling required to calculate TI and SS is straightforward.
However, we use the Metropolis algorithm to sample these densities and thus
allow a certain degree of autocorrelation in the samples, making the analysis
more realistic in an evidence estimation context. The Markov chains have a lag
of around 18 on average. In addition, we consider independent samples to
assess MBB performance in the context of error estimation.

(a) $\beta$ values spread according to evenly spaced quantiles of a Beta(0.3,
1) distribution.

(b) $\beta$ values equally spaced between 0 and 1.

Figure 1: Log-marginal likelihood estimates as a function of the number of
temperatures $K$ for the Gaussian model. Error bars depict $\pm 1$ standard
error based on 1000 independent MCMC analyses.

###  IV.1 Evidence estimate

We consider the following model specifications: $v=0.01$ and $d=20$. This
yields a log-marginal likelihood value $-46.15$. The analysis is performed for
$n=1000$ and $K=4,8,16,32,64,128,256$. Strictly speaking, SS uses $K-1$
temperatures, since it does not require samples from the posterior. For the
arrangement of the $\beta$ values, we test two approaches: evenly spaced
values from 0 to 1, and values spread according to evenly spaced quantiles of
a $\text{Beta}(0.3,1)$ distribution. The MCMC analysis is replicated 1000
times (with different random seeds) in order to calculate the error associated
with the estimates. The same power posterior samples are used to estimate SS
and TI.

Figures 1(a) and 1(b) display the results. It becomes clear in both cases that
the SS algorithm requires less temperatures than TI to produce estimates
around the true value. When the $\beta$ values are calculated according to a
Uniform(0,1), Figure 1(b), the TI estimates are seriously biased for low
number of temperatures, whereas the SS estimates, even though biased too, are
closer to the true value. For equally spaced $\beta$ values and $K=4$ Figure
1(b), TI is more than 130 units away from the true value compared to the
around 25 units for $\beta$ values spread according to quantiles of the
Beta(0.3,1) distribution in Figure 1(a). This shows that TI is more sensitive
to the distribution of the temperatures as was similarly shown by Xie  _et
al._ (2011).

Both methods improve their performance when most of the computational effort
is allocated in sampling in power posterior distributions near the prior, that
is for high temperatures. This is the effect of the Beta(0.3, 1) distribution,
which allows that half of the $\beta$ values are less than 0.1. The results
for this case are displayed in Figure 1(a). Even though TI improves its
performance considerably, it can not outperform SS, which still needs fewer
step temperatures to produce estimates around the true value.

###  IV.2 Standard error estimate

Based on the case that the $\beta$ values follow a Beta(0.3, 1) distribution,
we study the performance of the MBB method for estimating the evidence error.
For this, we calculate the standard error from the 1000 independent evidence
estimates used in the previous analysis, call this
$\widehat{\mbox{s.e.}}_{\text{ind}}$ and compare it to the standard error
estimates calculated using MBB ($\widehat{\mbox{s.e.}}_{\text{MBB}}$ for
different block lengths, $\ell=1,10,30,50,100,200,300$, via their differences,
$\widehat{\mbox{s.e.}}_{\text{MBB}}-\widehat{\mbox{s.e.}}_{\text{ind}}$.

(a) The Markov chains contain a degree of autocorrelation.

(b) The samples in the Markov chain are completely independent.

Figure 2: Difference between the standard error calculated via MBB and the one
from independent evidence estimates. The legend shows the different block
lengths used in MBB.

The results are shown in Figure 2(a). The case $\ell=1$ is the original
bootstrap method, which is used frequently for power posterior methods, but
which ignores the dependence in the sampled values of the Markov chain. It is
obvious that in the simple bootstrap with block length $\ell=1$, the standard
error is severely underestimated. On the other hand, the standard error
estimates improved significantly using the MBB with larger block lengths, but
still some underestimate the standard error. However, this example is an
extreme case of highly correlated Markov chains.

We have also performed the analysis in the ideal case that the samples in the
Markov chains are completely independent. The result are displayed in Figure
2(b). In this case, the standard bootstrap method, that is $\ell=1$, is
sufficient to estimate the standard error reasonably well. Large block lengths
cause, in general, a slight underestimation but only in the case of a low
number of temperatures. As the number of temperatures increases, the estimates
are located around the empirical error estimates, i.e., around zero, and less
dispersed.

We caution against the use of the theoretical standard error estimate of NS in
Equation (7) when the Metropolis-Hastings algorithm is used to generate the
samples rather than sampling independently, as the validity of this
theoretical standard error estimate is based on the independence assumption.
To this end, we include a comparison of this theoretical NS standard error
estimate with the empirical standard error obtained from 100 independent runs
in Table I. We observe a decrease in bias with increasing number of MCMC
steps. However, the NS standard error estimates are far too small and thus
underestimate the uncertainty even for a large number of MCMC steps of 5000.
This is a well known shortcoming of the NS standard error estimate, e.g. a
more detailed examination of this issue can also be found in Figure 4 of
Veitch and Vecchio (2010).

Table 1:  Nested sampling results based on 100 independent NS runs. $N$ is the number of live points, “Steps” the number of MCMC steps used to generate the points at each iteration, $Ave(\widehat{s.e.}_{\text{NS}})$ the average of the theoretical standard error estimate defined in (7), SD($\widehat{s.e.}_{\text{NS}}$) the standard deviation of the theoretical standard error estimates, $\widehat{s.e.}_{\text{NS,ind}}$ the standard error estimate based on the independent marginal likelihood estimates, and “Bias” the difference between the true value and the mean of the NS marginal likelihood estimates. $N$ | Steps | $Ave(\widehat{s.e.}_{\text{NS}})$ | SD($\widehat{s.e.}_{\text{NS}}$) | $\widehat{s.e.}_{\text{NS,ind}}$ | Bias
---|---|---|---|---|---
10 | 10 | 1.9922 | 0.0921 | 3.6244 | 4.7724
10 | 100 | 1.8918 | 0.0562 | 2.3078 | -0.3310
10 | 1000 | 1.8959 | 0.0542 | 2.1083 | -0.1580
10 | 5000 | 1.8939 | 0.0562 | 2.3454 | -0.3233

##  V Application with Simulated LIGO-Virgo Data

We apply the SS algorithm to an example analysis of a simulated binary black
hole coalescence signal in the Advanced LIGO Aasi  _et al._ (2015) and
Advanced Virgo Acernese  _et al._ (2015) gravitational wave detectors,
operating at design sensitivity. The data contained 4 s of simulated Gaussian
noise, generated using the design sensitivity curves of two Advanced LIGO
detectors (Hanford, Livingston) and the Advanced Virgo detector, plus the GW
signal. The simulated black hole binary had component masses 25
$\text{M}_{\odot}$and 13 $\text{M}_{\odot}$, and lay at a luminosity distance
of 614 Mpc, with a combined signal-to-noise ratio of 17.9 in the three-
detector network. The analysis was performed in the frequency range 40–512 Hz
using the IMRPhenomPv2 waveform approximant Hannam  _et al._ (2014). The
system’s total angular momentum was inclined at 95∘ to the line-of-sight to
the binary, and the primary and secondary black holes had dimensionless spin
magnitudes of 0.67 and 0.12, tilted at $45^{\circ}$ and $90^{\circ}$ to the
orbital angular momentum. This configuration produces a precession of the
orbital plane which results in a waveform that is not well approximated by a
non-spinning signal. The analysis was performed using the 15-dimensional
parameterised model for a quasi-circular black hole binary commonly used in
LIGO-Virgo analyses (e.g. Abbott  _et al._ (2016b, d)), implemented in the
LALInference package Veitch  _et al._ (2015).

We estimate the marginal likelihood via NS, TI and SS. For NS, we performed 32
runs with 2000 live points each. For TI and SS, we considered 31 temperatures,
evenly spaced on a logarithmic scale, with 4700 samples from each. From these
simulations, we ran TI and SS to estimate the evidence for $K=7,11,16$, and
$31$. To compute the standard error of the evidence estimates, we applied the
MBB method for different block lengths and took the one that yielded the
maximum standard deviation as a conservative way of estimation. The results
are displayed in Table 2 and visualized in Figure 3. The evidence estimates of
SS and NS are closer than of SS and TI. In the light of the performance of TI
in the simulation study, it seems that TI would have needed more temperatures
to achieve an evidence estimate consistent with both SS and NS. The standard
error of the NS evidence estimate is quite large, especially given that it was
calculated for a large number of live points.

Table 2:  Evidence estimates and corresponding standard errors, and the Bayes factor from the NS, TI and SS methods for different number of temperatures $K$. Method | $K$ | $\log z$ | SD | $\log B$
---|---|---|---|---
NS | - | -5730.82 | 0.36 | 103.33
TI | 7 | -5732.79 | 0.40 | 101.36
| 11 | -5732.23 | 0.32 | 101.92
| 16 | -5731.80 | 0.32 | 102.35
| 31 | -5731.52 | 0.27 | 102.63
SS | 7 | -5729.48 | 0.32 | 104.67
| 11 | -5730.14 | 0.14 | 104.01
| 16 | -5730.10 | 0.13 | 104.05
| 31 | -5730.15 | 0.13 | 104.00
Figure 3: Evidence estimates $\pm$ 2 standard errors. Subscripts in TI and SS
stand for the number of temperatures.

##  VI Discussion

SS is a method to estimate the marginal likelihood which has enjoyed great
popularity in phylogenetics where it has been shown to work well. It requires
less computational effort than TI to yield an accurate estimate of the
evidence. In a simulation study with a simple Gaussian model, we have shown
that it is less sensitive to the dispersal of the inverse temperature $\beta$
values and achieves a higher accuracy with a smaller number of power posterior
distributions. To the best of our knowledge, it has not been applied for
calculating the evidence of gravitational wave models yet. Its implementation
in this context should be straightforward since its main complexity lies with
sampling from the power posterior, like TI. However, this can be done by using
the parallel tempering method, which has been widely implemented in GW
software packages such as LALInference.

The performance of SS depends on its specifications, such as the number of
MCMC steps $n$ in each parallel tempering chain, the number of temperatures
$K$ and the distribution of the inverse temperature $\beta$ values. In
addition, it depends on how different the prior and the posterior are. To
mitigate the dependence on the prior distribution, we aim to explore a recent
extension of SS known as generalized steppingstone sampling (GSS; Fan  _et
al._ , 2011). This method makes use of a reference distribution which aims to
shorten the distance between the prior and the posterior. Even though it
requires posterior samples to construct the reference distribution, it could
be more accurate than its simple version and require less steps to yield the
same accuracy. For this, the reference distribution needs to be a reasonable
approximation of the posterior, otherwise it can dramatically fail (Maturana
Russel  _et al._ , 2018).

One of the drawbacks of power posterior methods is the lack of a direct
formula for the standard error of the evidence estimate. In practice, the
methods are run multiple times in order to obtain an empirical standard error
estimate. This brute force approach will prove too computationally expensive
in most practial applications. Alternatively, the standard bootstrap has been
applied. It is computationally much cheaper than the brute force approach, but
it does not take the dependencies within and between the Markov chains into
account. In this paper, we have proposed a moving block bootstrap method. This
approach has the ability to allow for potential autocorrelation within the
chains and cross-correlation between chains. We showed in Example IV.2 of our
simulation study that the standard bootstrap severely underestimates the
standard error in the presence of autocorrelation in Markov chains but that
the moving block bootstrap significantly improves the standard error estimates
of the evidence.

###### Acknowledgements.

We thank Claudia Kirch for helpful discussions on the block-bootstrap for time
series. We also thank the New Zealand eScience Infrastructure (NeSI) for their
high performance computing facilities, and the Centre for eResearch at the
University of Auckland for their technical support. The Observatoire de la
Côte d’Azur also provided support for this research. PM’s and RM’s work is
supported by Grant 3714568 from the University of Auckland Faculty Research
Development Fund and the DFG Grant KI 1443/3-1. JV is supported by STFC grant
ST/K005014/1. NC’s work is supported by NSF grants PHY-1806990 and
PHY-1505373. This paper has been given LIGO Document Number P1800299. All
analysis was conducted in R, an open-source statistical software available on
CRAN (cran.r-project.org) and LALInference.

## References

  * Saha and Williams (1994) P. Saha and T. B. Williams, The Astronomical Journal 107, 1295 (1994).
  * Christensen and Meyer (1998) N. Christensen and R. Meyer, Phys. Rev. D 58, 082001 (1998).
  * Christensen  _et al._ (2001) N. Christensen, R. Meyer, L. Knox, and B. Luey, Classical and Quantum Gravity 18, 2677 (2001).
  * Knox  _et al._ (2001) L. Knox, N. Christensen, and C. Skordis, The Astrophysical Journal Letters 563, L95 (2001).
  * Sharma (2017) S. Sharma, Annual Review of Astronomy and Astrophysics 55, 213 (2017).
  * Hinshaw  _et al._ (2013) G. Hinshaw _et al._ , The Astrophysical Journal Supplement Series 208, 19 (2013).
  * Planck Collaboration  _et al._ (2014) Planck Collaboration, Ade, P. A. R., _et al._ , Astronomy & Astrophysics 571, A16 (2014).
  * Planck Collaboration  _et al._ (2016) Planck Collaboration, Ade, P. A. R., _et al._ , Astronomy & Astrophysics 594, A13 (2016).
  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a).
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 241102 (2016b).
  * (11) R. Meyer and N. Christensen, Significance 13, 20.
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017a).
  * Zhu  _et al._ (2017) J. Zhu, J. Chen, W. Hu, and B. Zhang, National Science Review 4, 627 (2017).
  * Goldstein  _et al._ (2017) A. Goldstein, P. Veres, E. Burns, M. S. Briggs, R. Hamburg, D. Kocevski, C. A. Wilson-Hodge, R. D. Preece, S. Poolakkil, O. J. Roberts, C. M. Hui, V. Connaughton, J. Racusin, A. von Kienlin, T. D. Canton, N. Christensen, T. Littenberg, K. Siellez, L. Blackburn, J. Broida, E. Bissaldi, W. H. Cleveland, M. H. Gibby, M. M. Giles, R. M. Kippen, S. McBreen, J. McEnery, C. A. Meegan, W. S. Paciesas, and M. Stanbro, The Astrophysical Journal Letters 848, L14 (2017).
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Classical and Quantum Gravity 32, 074001 (2015).
  * Acernese  _et al._ (2015) F. Acernese _et al._ , Classical and Quantum Gravity 32, 024001 (2015).
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ , The Astrophysical Journal Letters 848, L12 (2017b).
  * Singer and Price (2016) L. P. Singer and L. R. Price, Phys. Rev. D 93, 024013 (2016).
  * Veitch  _et al._ (2015) J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, W. Del Pozzo, F. Feroz, J. Gair, C.-J. Haster, V. Kalogera, T. Littenberg, I. Mandel, R. O’Shaughnessy, M. Pitkin, C. Rodriguez, C. Röver, T. Sidery, R. Smith, M. Van Der Sluys, A. Vecchio, W. Vousden, and L. Wade, Phys. Rev. D 91, 042003 (2015).
  * Abbott  _et al._ (2016c) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 6, 041014 (2016c).
  * Abbott  _et al._ (2016d) B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. X 6, 041015 (2016d), arXiv:1606.04856 [gr-qc] .
  * Abbott  _et al._ (2018a) B. P. Abbott _et al._ (Virgo, LIGO Scientific), (2018a), arXiv:1805.11579 [gr-qc] .
  * Abbott  _et al._ (2016e) B. P. Abbott _et al._ (LIGO Scientific and Virgo Collaborations), Phys. Rev. Lett. 116, 221101 (2016e).
  * Abbott  _et al._ (2018b) B. P. Abbott _et al._ (Virgo, LIGO Scientific), (2018b), arXiv:1805.11581 [gr-qc] .
  * Weinberg (2018) N. N. Weinberg (Virgo, LIGO Scientific), (2018), arXiv:1808.08676 [astro-ph.HE] .
  * Smith and Thrane (2018) R. Smith and E. Thrane, Phys. Rev. X 8, 021019 (2018).
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 141101 (2017c).
  * Callister  _et al._ (2017) T. Callister, A. S. Biscoveanu, N. Christensen, M. Isi, A. Matas, O. Minazzoli, T. Regimbau, M. Sakellariadou, J. Tasson, and E. Thrane, Phys. Rev. X 7, 041058 (2017).
  * Abbott  _et al._ (2018c) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 120, 201102 (2018c).
  * Newton and Raftery (1994) M. A. Newton and A. E. Raftery, J. Roy. Statist. Soc. Ser. B 56, 3 (1994).
  * Xie  _et al._ (2011) W. Xie, P. O. Lewis, Y. Fan, L. Kuo, and M.-H. Chen, Syst. Biol. 60, 150 (2011).
  * Lartillot and Philippe (2006) N. Lartillot and H. Philippe, Syst. Biol. 55, 195 (2006).
  * Friel and Pettitt (2008) N. Friel and A. N. Pettitt, J. Roy. Stat. Soc. B 70, 589 (2008).
  * Neal (2001) R. M. Neal, Stat. Comput. 11, 125 (2001).
  * Skilling (2006) J. Skilling, Bayesian Analysis 1, 833 (2006).
  * Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys. Rev. D 81, 062003 (2010).
  * Brewer and Donovan (2015) B. J. Brewer and C. P. Donovan, Mon. Not. R. Astron. Soc. 448, 3206 (2015).
  * Feroz  _et al._ (2009) F. Feroz, M. Hobson, and M. Bridges, Monthly Notices of the Royal Astronomical Society 398, 1601 (2009), cited By 723.
  * Henderson  _et al._ (2017) R. Henderson, P. Goggans, and L. Cao, Digital Signal Processing: A Review Journal 70, 84 (2017), cited By 1.
  * Maturana Russel  _et al._ (2018) P. Maturana Russel, B. J. Brewer, S. Klaere, and R. R. Bouckaert, Syst. Biol. , syy050 (2018).
  * Maturana Russel (2018) P. Maturana Russel, in  _Bayesian Inference and Maximum Entropy Methods in Science and Engineering_ , edited by A. Polpo, J. Stern, F. Louzada, R. Izbicki, and H. Takada (Springer International Publishing, Cham, 2018) pp. 211–219.
  * Green (1995) P. J. Green, Biometrika 82, 711 (1995).
  * Umstätter  _et al._ (2005) R. Umstätter, N. Christensen, M. Hendry, R. Meyer, V. Simha, J. Veitch, S. Vigeland, and G. Woan, Phys. Rev. D 72, 022001 (2005).
  * Cornish and Littenberg (2015) N. J. Cornish and T. B. Littenberg, Class. Quant. Grav. 32, 135012 (2015).
  * Gelman and Meng (1998) A. Gelman and X. Meng, Statistical Science 13, 163 (1998), cited By 437\.
  * Neal (1993) R. M. Neal, (1993).
  * Maturana Russel (2017) P. Maturana Russel,  _Bayesian inference in phylogenetics using Nested Sampling_ , Ph.D. thesis, The University of Auckland (2017).
  * Efron (1979) B. Efron, Ann. Statist. 7, 1 (1979).
  * Kunsch (1989) H. R. Kunsch, Ann. Statist. 17, 1217 (1989).
  * Lahiri (2003) S. N. Lahiri,  _Resampling Methods for Dependent Data_ , Springer series in statistics (Springer New York, New York, NY, 2003).
  * Hannam  _et al._ (2014) M. Hannam _et al._ , Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Fan  _et al._ (2011) Y. Fan, R. Wu, M.-H. Chen, L. Kuo, and P. O. Lewis, Mol. Biol. Evol. 28, 523 (2011).