Marginal Likelihood from the Gibbs Output

Summary

Chib (1995) derives a stable, automatic method for computing the marginal likelihood — the key ingredient for Bayes factors — directly from the output of a Gibbs sampler. The approach exploits the Basic Marginal Likelihood Identity (BMI): the marginal likelihood equals the likelihood times the prior divided by the posterior ordinate, which holds for any parameter value. Because the posterior ordinate at a high-density point can be estimated by Rao-Blackwellizing over the Gibbs draws, the entire calculation is a by-product of the Markov chain Monte Carlo (MCMC) run already performed for estimation. Applications to binary probit variable selection, Gaussian mixture models, and a Markov switching model demonstrate accuracy and negligible numerical standard errors.

Key Claims

• BMI identity: $m(\mathbf{y}) = f(\mathbf{y}|\boldsymbol{\theta}) \pi(\boldsymbol{\theta}) / \pi(\boldsymbol{\theta}|\mathbf{y})$ holds for any $\boldsymbol{\theta}$; evaluated in log scale at a high-density point $\boldsymbol{\theta}^*$, it gives $\ln\hat{m}(\mathbf{y}) = \ln f(\mathbf{y}|\boldsymbol{\theta}^*) + \ln\pi(\boldsymbol{\theta}^*) - \ln\hat{\pi}(\boldsymbol{\theta}^*|\mathbf{y})$.
• The posterior ordinate $\hat{\pi}(\boldsymbol{\theta}^*|\mathbf{y})$ is estimated by averaging the full conditional density at $\boldsymbol{\theta}^*$ over all Gibbs draws — a simulation-consistent Rao-Blackwell estimate.
• For $B$-block models, $\pi(\boldsymbol{\theta}^*|\mathbf{y})$ is factorized sequentially as $\prod_{r=1}^B \pi(\boldsymbol{\theta}_r^*|\mathbf{y}, \boldsymbol{\theta}_1^*, \ldots, \boldsymbol{\theta}_{r-1}^*)$; each reduced conditional ordinate requires a short additional Gibbs run with earlier blocks fixed at $\boldsymbol{\theta}^*$.
• Numerical standard errors (NSEs) are derived from the spectral density matrix of the Gibbs output using the delta method and a Newey-West (1987) lag-window estimator.
• Harmonic mean estimator (Newton-Raftery 1994) is simulation-consistent but unstable (inverse-likelihood lacks finite variance); this method is both consistent and stable.
• Binary probit application (53 patients, 9 models): marginal likelihood precisely estimated (NSE ≤ 0.024); model $M_8$ ($C + \log x_2 + x_3 + x_4$) wins, consistent with Collett (1991).
• Galaxy velocity mixture application: three-component unrestricted variance model ($\log m = -224.138$) preferred over restricted two-component ($-240.464$).
• Markov switching gross national product (GNP) application: switching model ($\log m = -229.496$) preferred over a first-order autoregressive AR(1) model ($-231.94$).

Concepts Introduced or Extended

• Marginal Data Density — the BMI identity and B-block reduced conditional algorithm are the canonical Chib (1995) method
• Gibbs Sampler — Rao-Blackwell posterior ordinate estimator and numerical standard error derivation
• Binary Probit — variable selection via Bayes factors across 9 models
• Mixture of Normals — determining number of components via marginal-data-density (MDD) comparison

Entities Mentioned

• Siddhartha Chib

Quotes

"The approach is fully automatic and stable, requiring no inputs beyond the draws from the simulation. Thus draws from the prior, or additional maximizations, or importance sampling functions, or any other tuning function, are not required."

"It was shown that the numerical standard error of the estimate can be derived from the posterior sample and the calculations are exhibited in problems dealing with probit regression and finite-mixture models."

My Take

This is the canonical reference for MCMC-based marginal likelihood computation. The BMI identity is a simple algebraic rearrangement of Bayes' theorem, but exploiting it for marginal likelihood estimation required the insight that Gibbs output contains enough information to estimate the posterior ordinate accurately. The key limitation — that the full conditionals must have closed-form expressions to allow Rao-Blackwellization — was later addressed by Chib and Jeliazkov (2001) for Metropolis-Hastings (M-H) samplers. The three B-block extension is straightforward in principle but adds B-1 reduced runs in practice; for large $B$ or expensive likelihoods, the Gelfand-Dey (1994) modified importance sampler or the Laplace approximation may be preferred (see Amisano-Federico (2004)). The NSE framework is a valuable addition: it makes the accuracy of the marginal likelihood estimate itself quantifiable.