Bayesian model comparison with un-normalised likelihoods

arXiv (Cornell University), 2009

Computation of the marginal likelihood from a simulated posterior distribution is central to Bayesian model selection but is computationally difficult. The often-used harmonic mean approximation uses the posterior directly but is unstably sensitive to samples with anomalously small values of the likelihood. The Laplace approximation is stable but makes strong, and often inappropriate, assumptions about the shape of the posterior distribution. It is useful, but not general. We need algorithms that apply to general distributions, like the harmonic mean approximation, but do not suffer from convergence and instability issues. Here, I argue that the marginal likelihood can be reliably computed from a posterior sample by careful attention to the numerics of the probability integral. Posing the expression for the marginal likelihood as a Lebesgue integral, we may convert the harmonic mean approximation from a sample statistic to a quadrature rule. As a quadrature, the harmonic mean approximation suffers from enormous truncation error as consequence . This error is a direct consequence of poor coverage of the sample space; the posterior sample required for accurate computation of the marginal likelihood is much larger than that required to characterize the posterior distribution when using the harmonic mean approximation. In addition, I demonstrate that the integral expression for the harmonic-mean approximation converges slowly at best for high-dimensional problems with uninformative prior distributions. These observations lead to two computationally-modest families of quadrature algorithms that use the full generality sample posterior but without the instability. The first algorithm automatically eliminates the part of the sample that contributes large truncation error. The second algorithm uses the posterior sample to assign probability to a partition of the sample space and performs the marginal likelihood integral directly. This eliminates convergence issues. The first algorithm is analogous to standard quadrature but can only be applied for convergent problems. The second is a hybrid of cubature: it uses the posterior to discover and tessellate the subset of that sample space was explored and uses quantiles to compute a representative field value. Qualitatively, the first algorithm improves the harmonic mean approximation using numerical analysis, and the second algorithm is an adaptive version of the Laplace approximation. Neither algorithm makes strong assumptions about the shape of the posterior distribution and neither is sensitive to outliers. Based on numerical tests, we recommend a combined application of both algorithms as consistency check to achieve a reliable estimate of the marginal likelihood from a simulated posterior distribution.

Statistics and Computing, 2015

Recent decades have seen enormous improvements in computational inference for statistical models; there have been competitive continual enhancements in a wide range of computational tools. In Bayesian inference, first and foremost, MCMC techniques have continued to evolve, moving from random walk proposals to Langevin drift, to Hamiltonian Monte Carlo, and so on, with both theoretical and algorithmic innovations opening new opportunities to practitioners. However, this impressive evolution in capacity is confronted by an even steeper increase in the complexity of the datasets to be addressed. The difficulties of modelling and then handling ever more complex datasets most likely call for a new type of tool for computational inference that dramatically reduces the dimension and size of the raw data while capturing its essential aspects. Approximate models and algorithms may thus be at the core of the next computational revolution.

Water Resources Research, 2017

What is the “best” model? The answer to this question lies in part in the eyes of the beholder, nevertheless a good model must blend rigorous theory with redeeming qualities such as parsimony and quality of fit. Model selection is used to make inferences, via weighted averaging, from a set of K candidate models, , and help identify which model is most supported by the observed data, . Here, we introduce a new and robust estimator of the model evidence, , which acts as normalizing constant in the denominator of Bayes’ theorem and provides a single quantitative measure of relative support for each hypothesis that integrates model accuracy, uncertainty, and complexity. However, is analytically intractable for most practical modeling problems. Our method, coined GAussian Mixture importancE (GAME) sampling, uses bridge sampling of a mixture distribution fitted to samples of the posterior model parameter distribution derived from MCMC simulation. We benchmark the accuracy and reliability ...

Bayesian Analysis, 2012

Computation of the marginal likelihood from a simulated posterior distribution is central to Bayesian model selection but is computationally difficult. The often-used harmonic mean approximation uses the posterior directly but is unstably sensitive to samples with anomalously small values of the likelihood. The Laplace approximation is stable but makes strong, and often inappropriate, assumptions about the shape of the posterior distribution. It is useful, but not general. We need algorithms that apply to general distributions, like the harmonic mean approximation, but do not suffer from convergence and instability issues. Here, I argue that the marginal likelihood can be reliably computed from a posterior sample by careful attention to the numerics of the probability integral. Posing the expression for the marginal likelihood as a Lebesgue integral, we may convert the harmonic mean approximation from a sample statistic to a quadrature rule. As a quadrature, the harmonic mean approximation suffers from enormous truncation error as consequence. This error is a direct consequence of poor coverage of the sample space; the posterior sample required for accurate computation of the marginal likelihood is much larger than that required to characterize the posterior distribution when using the harmonic mean approximation. In addition, I demonstrate that the integral expression for the harmonic-mean approximation converges slowly at best for high-dimensional problems with uninformative prior distributions. These observations lead to two computationally-modest families of quadrature algorithms that use the full generality sample posterior but without the instability. The first algorithm automatically eliminates the part of the sample that contributes large truncation error. The second algorithm uses the posterior sample to assign probability to a partition of the sample space and performs the marginal likelihood integral directly. This eliminates convergence issues. The first algorithm is analogous to standard quadrature but can only be applied for convergent problems. The second is a hybrid of cubature: it uses the posterior to discover and tessellate the subset of that sample space was explored and uses quantiles to compute a representative field value. Qualitatively, the first algorithm improves the harmonic mean approximation using numerical analysis, and the second algorithm is an adaptive version of the Laplace approximation. Neither algorithm makes strong assumptions about the shape of the posterior distribution and neither is sensitive to outliers. Based on numerical tests, we recommend a combined application of both algorithms as consistency check to achieve a reliable estimate of the marginal likelihood from a simulated posterior distribution.

Data Mining and Knowledge Discovery, 2003

Markov chain Monte Carlo (MCMC) techniques revolutionized statistical practice in the 1990s by providing an essential toolkit for making the rigor and flexibility of Bayesian analysis computationally practical. At the same time the increasing prevalence of massive datasets and the expansion of the field of data mining has created the need for statistically sound methods that scale to these large problems. Except for the most trivial examples, current MCMC methods require a complete scan of the dataset for each iteration eliminating their candidacy as feasible data mining techniques. In this article we present a method for making Bayesian analysis of massive datasets computationally feasible. The algorithm simulates from a posterior distribution that conditions on a smaller, more manageable portion of the dataset. The remainder of the dataset may be incorporated by reweighting the initial draws using importance sampling. Computation of the importance weights requires a single scan of the remaining observations. While importance sampling increases efficiency in data access, it comes at the expense of estimation efficiency. A simple modification, based on the “rejuvenation” step used in particle filters for dynamic systems models, sidesteps the loss of efficiency with only a slight increase in the number of data accesses. To show proof-of-concept, we demonstrate the method on two examples. The first is a mixture of transition models that has been used to model web traffic and robotics. For this example we show that estimation efficiency is not affected while offering a 99% reduction in data accesses. The second example applies the method to Bayesian logistic regression and yields a 98% reduction in data accesses.

Bayesian Analysis, 2011

Statistical methods of inference typically require the likelihood function to be computable in a reasonable amount of time. The class of "likelihood-free" methods termed Approximate Bayesian Computation (ABC) is able to eliminate this requirement, replacing the evaluation of the likelihood with simulation from it. Likelihood-free methods have gained in efficiency and popularity in the past few years, following their integration with Markov Chain Monte Carlo (MCMC) and Sequential Monte Carlo (SMC) in order to better explore the parameter space. They have been applied primarily to estimating the parameters of a given model, but can also be used to compare models. Here we present novel likelihood-free approaches to model comparison, based upon the independent estimation of the evidence of each model under study. Key advantages of these approaches over previous techniques are that they allow the exploitation of MCMC or SMC algorithms for exploring the parameter space, and that they do not require a sampler able to mix between models. We validate the proposed methods using a simple exponential family problem before providing a realistic problem from human population genetics: the comparison of different demographic models based upon genetic data from the Y chromosome.

2007

This paper proposes a new versatile strategy for sampling posterior distributions for analytically intractable models. Building such samplers using Markov Chain Monte Carlo methodology usually leads to algorithms which are rather expensive from a computational point of view, hence hav- ing very few chances to be used for practical applications. The strategy we propose overcomes this drawback and is easy to use. A direct applica- tion of the proposed method is shown and discussed: maximum likelihood estimation of the parameters of a spatial point process. Resume Cet article propose une methodologie souple pour la simulation des lois a posteriori des modeles ayant des constantes de normalisation qui ne sont pas calculables analytiquement. Utiliser na¨õvement la philosophie Monte Carlo pour ce type de probleme, amenea des algorithmes couteux du point de vue du temps de calcul, donc peu utilisables en pratique. La strategie que nous proposonselimine cette difficulte et de plus elle est ...

SSRN Electronic Journal, 2000

This paper presents the R package MitISEM (mixture of t by importance sampling weighted expectation maximization) which provides an automatic and flexible two-stage method to approximate a non-elliptical target density kernel-typically a posterior density kernel-using an adaptive mixture of Student t densities as approximating density. In the first stage a mixture of Student t densities is fitted to the target using an expectation maximization algorithm where each step of the optimization procedure is weighted using importance sampling. In the second stage this mixture density is a candidate density for efficient and robust application of importance sampling or the Metropolis-Hastings (MH) method to estimate properties of the target distribution. The package enables Bayesian inference and prediction on model parameters and probabilities, in particular, for models where densities have multi-modal or other non-elliptical shapes like curved ridges. These shapes occur in research topics in several scientific fields. For instance, analysis of DNA data in bio-informatics, obtaining loans in the banking sector by heterogeneous groups in financial economics and analysis of education's effect on earned income in labor economics. The package MitISEM provides also an extended algorithm, 'sequential Mi-tISEM', which substantially decreases computation time when the target density has to be approximated for increasing data samples. This occurs when the posterior or predictive density is updated with new observations and/or when one computes model probabilities using predictive likelihoods. We illustrate the MitISEM algorithm using three canonical statistical and econometric models that are characterized by several types of non-elliptical posterior shapes and that describe well-known data patterns in econometrics and finance. We show that MH using the candidate density obtained by MitISEM outperforms, in terms of numerical efficiency, MH using a simpler candidate, as well as the Gibbs sampler. The MitISEM approach is also used for Bayesian model comparison using predictive likelihoods.

Journal of Statistical Planning and Inference, 2007

Many commonly used models in statistics can be formulated as (Bayesian) hierarchical Gaussian Markov random field (GMRF) models. These are characterised by assuming a (often large) GMRF as the second stage in the hierarchical structure and a few hyperparameters at the third stage. Markov chain Monte Carlo (MCMC) is the common approach for Bayesian inference in such models. The variance of the Monte Carlo estimates is O p (M −1/2 ) where M is the number of samples in the chain so, in order to obtain precise estimates of marginal densities, say, we need M to be very large.

ADVANCES IN NEURAL INFORMATION …, 2007