Computation of the marginal likelihood from MCMC samples

This is a recurring question (see this post, this post and this post), but I have a different spin.

Suppose I have a bunch of samples from a generic MCMC sampler. For each sample $\theta$, I know the value of the log likelihood $\log f(\textbf{x} | \theta)$ and of the log prior $\log f(\theta)$. If it helps, I also know the value of the log likelihood per data point, $\log f(x_i | \theta)$ (this information helps with certain methods, such as WAIC and PSIS-LOO).

I want to obtain a (crude) estimate of the marginal likelihood, just with the samples that I have, and possibly a few other function evaluations (but without rerunning an ad hoc MCMC).

First of all, let's clear the table. We all know that the harmonic estimator is the worst estimator ever. Let's move on. If you are doing Gibbs sampling with priors and posteriors in closed form, you can use Chib's method; but I am not sure how to generalize outside of those cases. There are also methods that require you to modify the sampling procedure (such as via tempered posteriors), but I am not interested in that here.

The approach I am thinking of consists of approximating the underlying distribution with a parametric (or nonparametric) shape $g(\theta)$, and then figuring out the normalization constant $Z$ as a 1-D optimization problem (i.e., the $Z$ that minimizes some error between $Z g(\theta)$ and $f(\textbf{x}|\theta) f(\theta)$, evaluated on the samples). In the simplest case, suppose that the posterior is roughly multivariate normal, I can fit $g(\theta)$ as a multivariate normal and get something similar to a Laplace approximation (I might want to use a few additional function evaluations to refine the position of the mode). However, I could use as $g(\theta)$ a more flexible family such as a variational mixture of multivariate $t$ distributions.

I appreciate that this method works only if $Z g(\theta)$ is a reasonable approximation to $f(\textbf{x}|\theta) f(\theta)$, but any reason or cautionary tale of why it would be very unwise to do it? Any reading that you would recommend?

The fully nonparametric approach uses some nonparametric family, such as a Gaussian process (GP), to approximate $\log f(\textbf{x}|\theta) + \log f(\theta)$ (or some other nonlinear transformation thereof, such as the square root), and Bayesian quadrature to implicitly integrate over the underlying target (see here and here). This seems to be an interesting alternative approach, but analogous in spirit (also, note that GPs would be unwieldy in my case).

• I think Chib, S. and Jeliazkov, I. 2001 "Marginal likelihood from the Metropolis--Hastings output" generalizes to normal MCMC outputs - would be interested to hear experiences with this approach. As for the GP - basically, this boils down to emulation of the posterior, which you could also consider for other problems. I guess the problem is that you are never sure about the quality of the approximation. What I also wonder is if an MCMC sample is ideal for a GP model, or whether you should invest more in the tails. – Florian Hartig Apr 29 '16 at 9:31
• (+1) Thanks for the reference, looks spot on -- I will check it out. I agree that all model-based approaches can be problematic (the good thing with Bayesian quadrature is that you get an estimate of uncertainty, although not sure how calibrated it is). For the moment my modest goal is to do something that is "better than a Laplace approximation". – lacerbi Apr 29 '16 at 15:48

The extension by Chib and Jeliazkov (2001) unfortunately gets quickly costly or highly variable, which is a reason why it is not much used outside Gibbs sampling cases.

While there are many ways and approaches to the normalisation constant $\mathfrak{Z}$ estimation problem (as illustrated by the quite diverse talks in the Estimating Constant workshop we ran last week at the University of Warwick, slides available there), some solutions do exploit directly the MCMC output.

1. As you mentioned, the harmonic mean estimator of Newton and Raftery (1994) is almost invariably poor for having an infinite variance. However, there are ways to avoid the infinite variance curse by using instead a finite support target in the harmonic mean identity $$\int \dfrac{\alpha(\theta)}{\pi(\theta)f(x|\theta)}\text{d}\pi(\theta|x)=\frac{1}{\mathfrak{Z}}$$ by picking $\alpha$ as the indicator of an HPD region for the posterior. This ensures finite variance by removing the tails in the harmonic mean. (Details are to be found in a paper I wrote with Darren Wraith and in a chapter about normalising constants written with Jean-Michel Marin.) In short, the method recycles the MCMC output $\theta_1,\ldots,\theta_M$ by identifying the $\beta$ (20% say) largest values of the target $\pi(\theta)f(x|\theta)$ and creating $\alpha$ as a uniform over the union of the balls centred at those largest density (HPD) simulations $\theta^0_i$ and with radius $\rho$, meaning the estimate of the normalising constant $\mathfrak{Z}$ is given by $$\hat{\mathfrak{Z}}^{-1}=\underbrace{\frac{1}{\beta M^2}\sum_{m=1}^M}_{\text{double sum over}\\\beta M\text{ ball centres }\theta_i^0\\\text{and M simulations } \theta_m} \underbrace{\mathbb{I}_{(0,\rho)}(\min_i||\theta_m-\theta^0_i||)\{\pi(\theta_m)f(x|\theta_m)\}^{-1}\big/\overbrace{\pi^{d/2}\rho^d\Gamma(d/2+1)^{-1}}^{\text{volume of ball with radius \rho}}}_{\dfrac{\beta M\alpha(\theta_m)}{\pi(\theta_m)f(x|\theta_m)}}$$ if $d$ is the dimension of $\theta$ (corrections apply for intersecting balls) and if $\rho$ is small enough for the balls to never intersect (meaning that at best only one indicator on the balls is different from zero). The explanation for the $\alpha M^2$ denominator is that this is a double sum of $\beta M^2$ terms: $$\frac{1}{\beta M}\sum_{i=1}^{\beta M} \underbrace{\frac{1}{M}\sum_{m=1}^M {\cal U}(\theta_i^0,\rho)(\theta_m)}_{\text{same as with \min}} \times \frac{1}{\pi(\theta_m)f(x|\theta_m)}$$ with each term in $\theta_m$ integrating to ${\mathfrak{Z}}^{-1}$.

2. Another approach is to turn the normalising constant $\mathfrak{Z}$ into a parameter. This sounds like a statistical heresy but the paper by Guttmann and Hyvärinen (2012) convinced me of the opposite. Without getting too much into details, the neat idea therein is to turn the observed log-likelihood $$\sum_{i=1}^n f(x_i|\theta) - n \log \int \exp f(x|\theta) \text{d}x$$ into a joint log-likelihood $$\sum_{i=1}^n[f(x_i|\theta)+\nu]-n\int\exp[f(x|\theta)+\nu]\text{d}x$$ which is the log-likelihood of a Poisson point process with intensity function $$\exp\{ f(x|\theta) + \nu +\log n\}$$ This is an alternative model in that the original likelihood does not appear as a marginal of the above. Only the modes coincide, with the conditional mode in ν providing the normalising constant. In practice, the above Poisson process likelihood is unavailable and Guttmann and Hyvärinen (2012) offer an approximation by means of a logistic regression. To connect even better with your question, Geyer's estimate is a MLE, hence solution to a maximisation problem.

3. A connected approach is Charlie Geyer's logistic regression approach. The fundamental notion is to add to the MCMC sample from $\pi(\theta|x)$ another sample from a known target, e.g., your best guess at $\pi(\theta|x)$, $g(\theta)$, and to run logistic regression on the index of the distribution behind the data (1 for $\pi(\theta|x)$ and 0 for $g(\theta)$). With the regressors being the values of both densities, normalised or not. This happens to be directly linked with Gelman and Meng (1997) bridge sampling, which also recycles samples from different targets. And later versions, like Meng's MLE.
4. A different approach that forces one to run a specific MCMC sampler is Skilling's nested sampling. While I [and others] have some reservations on the efficiency of the method, it is quite popular in astrostatistics and cosmology, with software available like multinest.
5. A last [potential if not always possible] solution is to exploit the Savage-Dickey representation of the Bayes factor in the case of an embedded null hypothesis. If the null writes as $H_0: \theta=\theta_0$ about a parameter of interest and if $\xi$ is the remaining [nuisance] part of the parameter of the model, assuming a prior of the form $\pi_1(\theta)\pi_2(\xi)$, the Bayes factor of $H_0$ against the alternative writes as $$\mathfrak{B}_{01}(x)=\dfrac{\pi^\theta(\theta_0|x)}{\pi_1(\theta_0)}$$ where $\pi^\theta(\theta_0|x)$ denotes the marginal posterior density of $\theta$ at the specific value $\theta_0$. In case the marginal density under the null $H_0: \theta=\theta_0$ $$m_0(x)=\int_\Xi f(x|\theta_0,\xi)\pi_2(\xi)\text{d}\xi$$ is available in closed form, one can derive the marginal density for the unconstrained model $$m_a(x)=\int_{\Theta\times\Xi} f(x|\theta,\xi)\pi_1(\theta)\pi_2(\xi)\text{d}\theta\text{d}\xi$$ from the Bayes factor. (This Savage-Dickey representation relies on specific versions of three different densities and so is fraught with danger, not even mentioning the computational challenge of producing the marginal posterior.)

[Here is a set of slides I wrote about estimating normalising constants for a NIPS workshop last December.]

• (+1) Incredibly rich answer, thank you. This will be useful to me and, I suppose, many other people. It will take me some time to have a look at the various approaches, and then I might come back with specific questions. – lacerbi Apr 30 '16 at 18:37
• Starting from point (1)... I read the relevant articles. The "corrected" harmonic mean estimator seems exactly what I was looking for. It's neat and easy to compute given a MCMC output. So... what's the catch? It doesn't look like the method is being widely used, judging from a quick search on Google Scholar. What are its limitations? (besides the need to identify the HPD regions, which I imagine might become an issue for very complicated posteriors in high dimension). I am definitely going to give it a try -- but I wonder if there is something I need to be wary of. – lacerbi May 1 '16 at 2:42
• I added a few more details: the issue in implementing the HPD uniform is to figure out a proper compact approximation for the HPD region. The convex hull of points with high posterior values is (NP?) hard to determine while balls centred at those points may intersect, which creates a secondary normalising constant problem. – Xi'an May 1 '16 at 7:40
• @Xi'an : very helpful, thanks! Can I ask: of all the mentioned approaches, what would currently be your recommendation if one looks for a general approach that tends to work out of the box (i.e. no tuning / checking required from the user)? I would be especially interested in the case of models with a low (< 50) number of parameters, non-normal posteriors, and strong correlations between parameters. – Florian Hartig May 4 '16 at 13:56
• @FlorianHartig: the fact that a generic software like BUGS does not return a generic estimate of $\mathfrak{Z}$ is sort of revealing the extent of the problem. The many solutions that one can find in the specialised literature have not produced a consensus estimate. Hence, my recommendation would be to opt for Geyer's logistic regression solution, which is somewhat insensitive to dimension. – Xi'an May 4 '16 at 14:26