marginal_likelihood#

cultivars.bayes.evidence.marginal_likelihood(result)[source]#

The log marginal likelihood of a fitted result, by whichever route it admits.

The marginal likelihood \(m(y) = \int f(y \mid \theta)\, \pi(\theta)\, d\theta\) is the probability the model assigned to the data before seeing them, and its logarithm is the one number a posterior can report in place of a log likelihood. This function is the dispatcher: it calls the result’s own marginal_likelihood() method, which knows the route its posterior admits, and hands back the record so that results of different kinds can be ranked by compare() or weighted by bayesian_model_average() without the caller knowing which estimator produced each number.

Three routes exist in the package. The conjugate BVAR integrates \((B, \Sigma)\) out in closed form, so its record is exact with a zero Monte Carlo error. The GibbsBVAR under a static prior uses Chib’s (1995) identity \(\log m(y) = \log f(y \mid \theta^*) + \log \pi(\theta^*) - \log \pi(\theta^* \mid y)\), with the posterior ordinate at a high-density point estimated from the Gibbs output and no extra runs. The particle-chain PerturbationDSGE posterior uses the modified harmonic mean of modified_harmonic_mean(), and inherits the noise of its particle-filter likelihood.

Parameters:

result (object) – A fitted result exposing a marginal_likelihood() method that returns a record: the three results above. Duck-typed, so a result from outside the package that returns the same record type is accepted.

Return type:

_MarginalLikelihoodSelection

Returns:

The MarginalLikelihoodSelection the result produced: log_value, its Monte Carlo standard error mcse (zero when exact), the method, the source label, and nobs, the number of observations the likelihood scored, which compare() uses to refuse comparisons across different effective samples.

Raises:

SpecificationError – If the result offers no marginal likelihood, or if its method returns something other than the record. The sampled results that do not offer one – the stochastic-volatility family, the Student-t and stochastic-volatility BVARs, the time-varying-parameter models, the Gibbs BVAR under an adaptive prior or with dummy observations – need reduced-run or particle-based estimators that are not implemented; compare those by out-of-sample predictive score through stacking() or a backtest.

See also

References

Chib, S. (1995). Marginal likelihood from the Gibbs output. Journal of the American Statistical Association, 90(432), 1313-1321.

Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773-795.

Example

The conjugate BVAR’s record is exact:

>>> import numpy as np
>>> from cultivars.multivariate.large_dim.bayesian import BVAR
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((120, 2))
>>> for t in range(1, 120):
...     y[t] = 0.7 * y[t - 1] + rng.standard_normal(2)
>>> record = marginal_likelihood(BVAR(y, order=1).fit(n_draws=100, seed=0))
>>> record.method, record.mcse, record.nobs
('analytic', 0.0, 119)
>>> record.source
'BVAR(1), niw(l1=0.2, l3=1, l4=100)'

A point estimate has no marginal likelihood and is refused by name:

>>> from cultivars.multivariate.reduced_form.vector_autoregression import VAR
>>> marginal_likelihood(VAR(y, order=1).fit())
Traceback (most recent call last):
    ...
cultivars.exceptions.SpecificationError: VARResult offers no marginal likelihood. ...