GibbsBVARResult.marginal_likelihood

GibbsBVARResult.marginal_likelihood()

Chib’s (1995) log marginal likelihood of the sample.

The two-block identity

\[\log m(y) = \log p(y \mid \beta^*, \Sigma^*) + \log p(\beta^*, \Sigma^*) - \log p(\Sigma^* \mid y) - \log p(\beta^* \mid \Sigma^*, y),\]

at the posterior mean \((\beta^*, \Sigma^*)\), with the covariance ordinate averaged over the retained coefficient draws and the coefficient ordinate evaluated exactly from its Gaussian conditional. The Monte-Carlo standard error of the average is on the record.

Return type:

_MarginalLikelihoodSelection

Returns:

The record; compare() on it ranks models fitted to the same sample.

Raises:

SpecificationError – Under an adaptive prior, whose extra Gibbs block would need reduced runs, or a dummy-observation prior, whose artificial rows the sampler treats as data.

Example

Two static priors on one sample, ranked by Bayes factor:

>>> import numpy as np
>>> from cultivars.bayes.priors import IndependentNormalWishartPrior, SpikeAndSlabPrior
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((200, 3))
>>> for t in range(1, 200):
...     y[t] = 0.6 * y[t - 1] + rng.standard_normal(3)
>>> loose = GibbsBVAR(y, order=1).fit(n_draws=300, n_burn=100, seed=0)
>>> tight = GibbsBVAR(
...     y, order=1, prior=IndependentNormalWishartPrior(tightness=0.05)
... ).fit(n_draws=300, n_burn=100, seed=0)
>>> evidence = loose.marginal_likelihood()
>>> evidence.method, bool(evidence.mcse < 0.1)
('chib', True)
>>> bool(evidence.log_value > tight.marginal_likelihood().log_value)
True
>>> adaptive = GibbsBVAR(y, order=1, prior=SpikeAndSlabPrior()).fit(
...     n_draws=60, n_burn=20, seed=0
... )
>>> adaptive.marginal_likelihood()
Traceback (most recent call last):
cultivars.exceptions.SpecificationError: Chib's estimator here covers the two-block ...