Stepping-stone sampling algorithm for calculating the evidence of gravitational wave models

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 byproduct. 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.

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. (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 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. (a) The Markov chains contain a degree of autocorrelation. (b) The samples in the Markov chain are completely independent.

Figure 3

Evidence estimates $\pm 2$ standard errors. Subscripts in NS stand for the number of live points $N$, and in TI and SS stand for the number of temperatures $K$.