Product-of-Posteriors Estimator

• The product-of-posteriors estimator is defined by constructing a proposal density as the product of marginal posteriors from independent parameter blocks for efficient marginal likelihood estimation.
• It leverages a single MCMC run by re-mixing block samples through permutation to generate approximate i.i.d. draws, simplifying importance-sampling implementation.
• Empirical results in regression, normal mixtures, and Poisson models demonstrate its accuracy, finite variance, and computational efficiency when using appropriate block design and marginal approximations.

Searching arXiv for the cited paper and closely related marginal-likelihood estimation work to ground the article in the literature. arxiv_search(query="(1311.0674) product of posteriors marginal likelihood Perrakis Ntzoufras Tsionas", max_results=5) The product-of-posteriors estimator is an importance-sampling estimator of the marginal likelihood in which the proposal density is formed by multiplying the marginal posterior distributions of parameter blocks rather than using the full joint posterior. In the formulation studied by Perrakis, Ntzoufras, and Tsionas, the method is designed for multi-block parameter vectors, does not require additional Markov Chain Monte Carlo sampling, and is not dependent on the type of MCMC scheme used to sample from the posterior. Its primary use is marginal-likelihood and Bayes-factor computation, with empirical illustrations in normal regression models, finite normal mixtures, and longitudinal Poisson models showing accurate marginal likelihood estimates (1311.0674).

1. Formal construction

Let $\mathcal D$ denote observed data and let $\theta$ be the full parameter vector of a model $M$. The marginal likelihood, or evidence, is

$$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$

where $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$.

Using the posterior identity $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$, one may write

$$p(\mathcal D) = \int \frac{p(\theta,\mathcal D)}{p(\theta\mid\mathcal D)} \,p(\theta\mid\mathcal D)\,\mathrm d\theta = \mathbb E_{p(\theta\mid\mathcal D)} \!\left[ \frac{p(\theta,\mathcal D)} {p(\theta\mid\mathcal D)} \right].$$

This is a zero-variance identity in principle, since $$p(\theta\mid\mathcal D)\propto p(\theta,\mathcal D)$$, but in practice the normalizing constant of the full posterior is not directly available.

The product-of-posteriors construction begins by partitioning the parameter vector into $B$ blocks,

$\theta = (\theta_1,\theta_2,\dots,\theta_B),$

and defining the importance density

$\theta$0

The proposal $\theta$1 retains the support of the posterior while imposing independence across blocks. The resulting importance-sampling identity is

$\theta$2

with weight

$\theta$3

Given an i.i.d. sample $\theta$4, the Monte Carlo estimator is

$\theta$5

which converges almost surely to $\theta$6 as $\theta$7.

2. Statistical properties and relation to alternative estimators

By construction,

$\theta$8

Its variance is

$\theta$9

provided that $M$0 (1311.0674). Finite variance is guaranteed when the tails of $M$1 are at least as heavy as those of the true joint posterior.

The method is situated among several standard marginal-likelihood estimators. The harmonic-mean estimator of Newton–Raftery is described as unbiased but typically having infinite variance. Bridge sampling in the sense of Meng–Wong is described as more stable but requiring tuning of a bridge function. Chib’s method is described as very accurate when reduced posterior ordinates are available, but as potentially requiring additional MCMC runs. Against this background, the product-of-posteriors estimator is characterized as requiring no extra MCMC, being easy to implement, and having finite variance in many practical settings, especially when parameters can be blocked so that block-wise posteriors are weakly dependent.

The regularity conditions stated for validity of the law of large numbers, the central limit theorem, and the finite-variance argument are explicit. The support of $M$2 must include the support of the joint posterior; the weight $M$3 must have finite second moment under $M$4; and mild integrability conditions, including continuity and boundedness in tails, are required. A common misunderstanding is to regard the product proposal as a surrogate for the full posterior itself. In fact, the construction is exact only at the level of marginal blocks; the dependence structure is deliberately removed.

3. Implementation from posterior simulation output

A notable feature of the method is that one can construct an approximate i.i.d. sample from

$M$5

using a single MCMC output from the full posterior. If $M$6 is a posterior chain, then for each block $M$7 the draws $M$8 are independently permuted, and the remixed parameter vector $M$9 is formed by taking block $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$0 from its corresponding permuted position. This produces draws that are approximately i.i.d. from the product of marginals. A simpler alternative is a single random permutation of each block, and if multiple MCMC chains are available they may be interleaved and permuted across chains.

The per-sample computational burden is determined by evaluation of the likelihood $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$1, the prior density $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$2, and the marginal posterior densities $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$3. Estimation of the latter is model- and sampler-dependent. When Gibbs sampling is used and full-conditional normalizing constants are known, Rao–Blackwell estimates are available. Otherwise, one may fit moment-matched normal or other simple parametric densities to each block, or use univariate kernel density estimates.

The implementation strategy is therefore modular. The original posterior sample is re-used, while the main approximation enters through the marginal-density estimates and the block factorization. This suggests that the estimator is particularly attractive in workflows where posterior simulation is already available and recomputation is costly.

4. Empirical behavior in model classes

The paper reports three principal empirical illustrations (1311.0674).

Model class	Configuration	Reported behavior
Normal linear regression	25 data points, 4 competing models	Estimates virtually indistinguishable from truth; MC error $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$4 on log-scale
Finite normal mixtures	Galaxy data, up to 4 components	With label-switching remedies, excellent agreement with Neal’s large-sample benchmarks; MC error $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$5–$$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$6
Longitudinal Poisson random-effects	Epileptic seizure counts, 59 subjects	3-block error $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$7; 4-block error $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$8; both compare well to Chib-based methods

In the normal linear regression example, the true log-marginal likelihood is available in closed form. The product-of-posteriors estimator, using either exact marginals or Rao–Blackwell marginals, gives estimates virtually indistinguishable from the truth, with Monte Carlo error approximately $$p(\mathcal D) = \int p(\mathcal D\mid \theta)\,\pi(\theta)\,\mathrm d\theta = \int p(\theta,\mathcal D)\,\mathrm d\theta,$$9 on the log scale. The paper further reports favorable comparison with Laplace–Metropolis, Chen’s estimator, Chib’s method, and optimal bridge sampling. It also states that the estimator captures changes under very diffuse $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$0-priors without re-running MCMC.

In finite normal mixtures, the main complication is label-switching. Because a $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$1-component mixture has $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$2 symmetric modes, a naive product estimator of the full posterior is biased low. Two remedies are reported: adding a $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$3 bias correction when posterior modes are well separated, and using random-permutation sampling so that all labelings are represented. Both remedies yield estimates in excellent agreement with Neal’s large-sample benchmarks, with Monte Carlo error between $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$4 and $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$5 on the log scale.

In the longitudinal Poisson random-effects example, the hierarchical model is organized into either three blocks or four blocks, with the latter including random effects explicitly. The 3-block estimator integrates out the random effects but requires inner importance sampling to approximate the high-dimensional integrated likelihood; its reported Monte Carlo error on the log scale is approximately $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$6. The 4-block estimator works directly with the full Poisson hierarchical likelihood but requires approximation of a 116-dimensional joint marginal; its reported error is approximately $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$7. Both versions compare well to the methods of Chib et al. (1998) and Chib and Jeliazkov (2001), for which errors of approximately $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$8 are reported.

5. Limitations, failure modes, and practical diagnostics

The principal limitation is dependence across blocks. If the block posteriors $$p(\theta,\mathcal D)=p(\mathcal D\mid\theta)\,\pi(\theta)$$9 are far from independent, then the product approximation $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$0 may place too little mass in the high-density regions of the true joint posterior, producing large variance. In such settings the proposal preserves support but can still be inefficient because it misrepresents posterior geometry.

A second limitation is the quality of the marginal approximations. Poor kernel fits or mis-specified normal approximations introduce bias in the denominator of the importance weights. This issue is structurally different from Monte Carlo variance: even with a large number of re-mixed draws, inaccurate marginal-density estimation can distort the estimator.

The paper recommends several practical checks. One is an empirical variance-check: double the sample size from $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$1 to $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$2 and examine whether the Monte Carlo error shrinks by roughly $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$3; failure of this pattern may indicate infinite variance of $$p(\theta\mid\mathcal D)=p(\theta,\mathcal D)/p(\mathcal D)$$4. Another is to compare the product-of-posteriors estimator with a pilot bridge-sampler or Chib estimate on a reduced run in order to gauge bias. For mixture models with label-switching, ex-post random-permutation sampling is recommended rather than fixed labeling or bias corrections. In hybrid settings where some full-conditionals are unavailable, Rao–Blackwell estimates may be combined with moment-matched normals or univariate kernels.

These limitations also clarify a common misconception. The estimator does not avoid the central difficulty of marginal-likelihood estimation; rather, it relocates the problem into the design of an importance density that is computationally convenient yet sufficiently close to the posterior to control variance.

6. Block design and methodological significance

Block partitioning is central to performance. The stated recommendation is to choose blocks so that within-block parameters are strongly dependent while between-block parameters are weakly dependent. Natural Gibbs blocks and orthogonal reparameterizations, as discussed in connection with Gilks and Roberts (1996) and Ghosh and Clyde (2011), are singled out as helpful for satisfying near-independence across blocks.

Methodologically, the estimator occupies an intermediate position between posterior-ordinate approaches and more general-purpose importance or bridge schemes. It re-uses one MCMC run, is applicable to latent-variable and multimodal settings, and can be adapted to multi-block parameterizations without requiring a specific MCMC architecture. This suggests that its main significance lies not in replacing all competing estimators, but in providing a computationally economical route to evidence estimation when posterior samples are already available and a sensible block factorization can be constructed.

The concluding characterization in the source is that the product-of-posteriors importance sampler is conceptually simple, re-uses one MCMC run, is often finite-variance and unbiased, is competitive with more elaborate methods, and is flexible for multi-block, latent-variable, or multimodal settings. With careful blocking and marginal fitting, it is presented as a robust tool for marginal-likelihood and Bayes-factor computation (1311.0674).