[Paper Review] Bayesian evidence: can we beat MultiNest using traditional MCMC methods?

Rutger van Haasteren|ArXiv.org|Nov 11, 2009

Pulsars and Gravitational Waves Research1 references3 citations

TL;DR

This paper presents a novel method to compute Bayesian evidence using standard MCMC chains without re-running simulations, leveraging a subset of uncorrelated samples to achieve high accuracy. It outperforms MultiNest in high-dimensional, peaked problems and reaches the theoretical Poissonian error limit when using uncorrelated chains.

ABSTRACT

Markov Chain Monte Carlo (MCMC) methods have revolutionised Bayesian data analysis over the years by making the direct computation of posterior probability densities feasible on modern workstations. However, the calculation of the prior predictive, the Bayesian evidence, has proved to be notoriously difficult with standard techniques. In this work a method is presented that lets one calculate the Bayesian evidence using nothing but the results from standard MCMC algorithms, like Metropolis-Hastings. This new method is compared to other methods like MultiNest, and greatly outperforms the latter in several cases. One of the toy problems considered in this work is the analysis of mock pulsar timing data, as encountered in pulsar timing array projects. This method is expected to be useful as well in other problems in astrophysics, cosmology and particle physics.

Motivation & Objective

• To develop a method for computing Bayesian evidence using only existing MCMC chains, avoiding the need for new sampling.
• To overcome the limitations of traditional methods like the harmonic mean estimator and parallel tempering, which suffer from inaccuracy or high computational cost.
• To assess whether standard MCMC methods can outperform advanced nested sampling algorithms like MultiNest in evidence computation.
• To evaluate the impact of MCMC chain correlation on evidence estimation accuracy and convergence.
• To provide a reliable, bootstrap-based error estimation for the computed Bayesian evidence.

Proposed method

• The method uses a subset of MCMC samples from the posterior distribution, specifically those within a region defined by a threshold on the posterior density.
• It estimates the Bayesian evidence as a weighted average over this subset, using the formula $ I = \frac{1}{N} \sum_{i \in F_t} \frac{\pi(\vec{\Theta}_i) L(\vec{\Theta}_i)}{p(\vec{\Theta}_i)} $, where $ F_t $ is the subset of samples above a threshold $ t $.
• The algorithm relies on storing normalized posterior values alongside parameter samples during standard MCMC runs.
• It applies a bootstrap resampling procedure to estimate the uncertainty of the evidence integral, ensuring reliable error estimation.
• The method is tested on both toy problems and a pulsar timing data model, comparing performance to MultiNest and parallel tempering.
• The approach is shown to reach the theoretical Poissonian error limit $ \sigma = I / \sqrt{N} $ when using uncorrelated chains, indicating optimal efficiency.

Experimental results

Research questions

• RQ1 Can Bayesian evidence be accurately computed using only standard MCMC chains without additional sampling?
• RQ2 How does the performance of this MCMC-based evidence estimator compare to MultiNest in high-dimensional, peaked posterior problems?
• RQ3 To what extent does MCMC chain correlation degrade the accuracy of Bayesian evidence estimation?
• RQ4 Can the method achieve the theoretical error limit of stochastic integration methods when using uncorrelated samples?
• RQ5 How does the choice of threshold $ c $ and sample size $ N $ affect the precision and stability of the evidence estimate?

Key findings

• The proposed method significantly outperforms MultiNest in terms of accuracy and efficiency for high-dimensional, peaked posterior distributions.
• When using uncorrelated MCMC samples, the method achieves the theoretical Poissonian error limit $ \sigma = I / \sqrt{N} $, indicating optimal performance.
• Correlated MCMC chains degrade the accuracy of the evidence estimate, requiring longer chains and increasing uncertainty, even with large $ N $.
• Increasing the number of samples $ N $ and the threshold $ c $ improves precision, but bootstrap error increases with higher $ c $, indicating reduced stability.
• The method is applicable to previously run MCMC chains, enabling evidence computation without re-running simulations.
• The algorithm is expected to be highly useful in astrophysics, cosmology, and particle physics for model selection tasks involving complex, high-dimensional models.