Bayesian inference for challenging scientific models

Abstract

Advances in technology and computation have led to ever more complicated scientific models of phenomena across a wide variety of fields. Many of these models present challenges for Bayesian inference, as a result of computationally intensive likelihoods, high-dimensional parameter spaces or large dataset sizes. In this thesis we show how we can apply developments in probabilistic machine learning and statistics to do inference with examples of these types of models.

As a demonstration of an applied inference problem involving a non-trivial likelihood computation, we show how a combination of optimisation and MCMC methods along with careful consideration of priors can be used to infer the parameters of an ODE model of the cardiac action potential.

We then consider the problem of pileup, a phenomenon that occurs in astronomy when using CCD detectors to observe bright sources. It complicates the fitting of even simple spectral models by introducing an observation model with a large number of continuous and discrete latent variables that scales with the size of the dataset. We develop an MCMC-based method that can work in the presence of pileup by explicitly marginalising out discrete variables and using adaptive HMC on the remaining continuous variables. We show with synthetic experiments that it allows us to fit spectral models in the presence of pileup without biasing the results. We also compare it to neural Simulation- Based Inference approaches, and find that they perform comparably to the MCMC-based approach whilst being able to scale to larger datasets.

As an example of a problem where we wish to do inference with extremely large datasets, we consider the Extreme Deconvolution method. The method fits a probability density to a dataset where each observation has Gaussian noise added with a known sample-specific covariance, originally intended for use with astronomical datasets. The existing fitting method is batch EM, which would not normally be applied to large datasets such as the Gaia catalog containing noisy observations of a billion stars. In this thesis we propose two minibatch variants of extreme deconvolution, based on an online variation of the EM algorithm, and direct gradient-based optimisation of the log-likelihood, both of which can run on GPUs. We demonstrate that these methods provide faster fitting, whilst being able to scale to much larger models for use with larger datasets.

We then extend the extreme deconvolution approach to work with non- Gaussian noise, and to use more flexible density estimators such as normalizing flows. Since both adjustments lead to an intractable likelihood, we resort to amortized variational inference in order to fit them. We show that for some datasets that flows can outperform Gaussian mixtures for extreme deconvolution, and that fitting with non-Gaussian noise is now possible.