Bayesian optimisation for Bayesian evidence (BOBE) -- a fast and efficient likelihood emulator for model selection

TL;DR

BOBE tackles the costly Bayesian evidence calculation by learning a Gaussian Process emulator of the full likelihood and guiding sampling through a log-evidence–aware acquisition (WIPStd). The method achieves up to ~$10^3$ fewer likelihood evaluations than traditional nested sampling while preserving accurate log-evidences and posterior inferences, with overhead decoupled from the likelihood time $t_L$. It demonstrates strong performance across toy, multimodal, and cosmological (ΛCDM) likelihoods and shows substantial speedups over PolyChord in several settings. This approach enables practical Bayesian model selection for expensive likelihoods, offering robust uncertainty quantification and scalable convergence diagnostics that integrate smoothly with cosmological pipelines such as Cobaya and JAX-based tooling.

Abstract

The formalism of Bayesian model selection provides a very elegant way of ranking different physical models in terms of how compatible they are with a given set of observed data. However, its practical application is often hampered by the challenge of having to compute the Bayesian evidence - a multi-dimensional integral over the product of likelihood and prior probability. This usually necessitates a large number of function calls to the likelihood, which may become prohibitive in case of "slow", costly to evaluate likelihoods. A possible solution to this problem lies in approximating the slow full likelihood by a fast emulated likelihood. In this paper, we introduce BOBE (Bayesian Optimisation for Bayesian Evidence), a method to construct a Gaussian Process Regression (GPR)-based emulator. BOBE utilises a Bayesian Optimisation algorithm designed specifically to (i) provide a realistic estimate of the emulator's uncertainty and its impact on the evidence calculation, and (ii) minimise the number of likelihood evaluations required in order to meet a given evidence accuracy goal. We apply it to a number of toy examples as well as actual cosmological likelihoods, and demonstrate that training the emulator to a sufficient accuracy takes a factor of $O(10^3)$ fewer direct likelihood evaluations than would be needed if one were to directly compute the evidence integral via nested sampling. BOBE's overhead is independent of the likelihood computation time $t_L$, making it particularly useful for "expensive" likelihoods with $t_L \gtrsim 1$~s. BOBE is written in Python, supports MPI parallelisation, takes advantage of automatic differentiation and just-in-time-compilation provided by JAX, can straightforwardly be implemented with cosmological data analysis frameworks such as Cobaya, and is available for download from https://github.com/Ameek94/BOBE.

Figures (13)

• Figure 1: Flowchart of the general BOBE framework for constructing a likelihood emulator and computing the log-evidence.
• Figure 2: Demonstration of the WIPStd acquisition function in action across three consecutive iterations. Top row: objective function (black solid line) and GP mean prediction (blue points) with evaluation locations (blue circles) and $\pm 2 \sigma(x)$ GP uncertainty region (blue bands). Bottom Row:$\mathrm{WIPStd}(x)$, approximated by $\widehat{\mathfrak{W}}(x)$ (light blue solid line) with evaluation locations (grey dashed lines) and location of next evaluation (red dashed line). Note the rapid decrease in $\mathrm{WIPStd}(x)$ as additional samples are added
• Figure 3: BOBE applied to a 2D Gaussian posterior. Top Left: Triangle plot comparing 1D and 2D posteriors contours obtained from the BOBE emulator (blue blobs/lines) and the true posterior (black lines). Locations of the training samples are indicated by red dots. Top Right: WIPStd convergence diagnostic (solid blue line) compared with the absolute difference between BOBE evidence estimate and the analytic evidence reference value (green dashed line), plotted as a function of the number of BO iterations. The red dashed line denotes the convergence criterion. Bottom Row: Grid evaluations of the true function (left), BOBE GP predictive mean (centre) and their difference (right). The run was initialised with 2 Sobol samples.
• Figure 4: Same as Figure \ref{['fig:gaussian']} for the Gaussian ring. This run was initialised with 8 Sobol samples.
• Figure 5: Same as Figure \ref{['fig:gaussian']} for the Himmelblau function. This run was initialised with 8 Sobol samples.