compare#

cultivars.bayes.evidence.compare(*records, names=None, prior_probabilities=None)[source]#

Rank models by marginal likelihood with posterior model probabilities.

One row per record, best first. Each row carries the log marginal likelihood, its Monte Carlo standard error, the log Bayes factor against the best model \(\log B = \log m_i(y) - \max_j \log m_j(y)\), the Kass-Raftery reading of \(2 \log B\) (from “not worth more than a bare mention” through “very strong”), the posterior model probability

\[p(M_i \mid y) = \frac{\pi_i\, m_i(y)}{\sum_j \pi_j\, m_j(y)},\]

and the method that produced the number. The notes warn when the records score different numbers of observations, since their differences are then not Bayes factors, and when a model sits within two Monte Carlo standard errors of the best, since the ranking is then not resolved by the estimates at hand.

Parameters:
Return type:

SummaryTable

Returns:

The _SummaryTable with columns model, log ML, mcse, log BF vs best, evidence, post. prob., and method; print it or call its to_pandas().

Raises:
  • SpecificationError – If fewer than two records are given, or a prior probability is negative or all are zero.

  • DimensionError – If names or prior_probabilities do not have one entry per record.

See also

References

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

Example

Two conjugate BVARs differing in prior tightness, ranked; the six-nat gap reads as very strong evidence for the looser prior:

>>> import numpy as np
>>> from cultivars.bayes.priors import NormalInverseWishartPrior
>>> 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)
>>> loose = marginal_likelihood(BVAR(y, order=1).fit(n_draws=100, seed=0))
>>> tight = BVAR(y, order=1, prior=NormalInverseWishartPrior(tightness=0.02))
>>> tight = marginal_likelihood(tight.fit(n_draws=100, seed=0))
>>> table = compare(loose, tight, names=("loose", "tight"))
>>> table.columns
('model', 'log ML', 'mcse', 'log BF vs best', 'evidence', 'post. prob.', 'method')
>>> table.rows[0]
('loose', '-357.408', '0.000', '0.000', 'best', '0.998', 'analytic')
>>> table.rows[1][0], table.rows[1][3], table.rows[1][4]
('tight', '-6.140', 'very strong')