bayesian_model_average#

cultivars.bayes.combination.bayesian_model_average(*results, records=None, prior_probabilities=None, names=None)[source]#

Weight models by posterior model probability.

For models \(M_1, \ldots, M_m\) fitted to the same sample \(y\), with marginal likelihoods \(p(y \mid M_i)\) and prior probabilities \(\pi_i\), the weight on model \(i\) is

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

computed on the log scale with a log-sum-exp normalization so that a spread of hundreds of nats between marginal likelihoods – routine for a VAR – does not underflow. The mixture predictive is then \(\sum_i w_i\, p(\tilde y \mid y, M_i)\), and the returned record draws from it by drawing a model in proportion to its weight and then one of that model’s predictive paths (Hoeting et al., 1999).

This is the coherent answer when one of the models is true, and as the sample grows the weights concentrate on it. When none is true the weights still concentrate, on whichever model is closest in Kullback-Leibler divergence, so the average collapses to a single model with extra steps; stacking() is the combination to prefer for forecasting under misspecification, and the record’s notes say so.

The marginal likelihoods must score the same observations. Two models of different autoregressive order condition on different presamples and therefore score samples of different length, and the ratio of their marginal likelihoods is then a ratio of densities of different data, not a Bayes factor. The function refuses that case; the remedy is to drop leading rows so every model conditions on the same presample and refit.

Parameters:

  • *results (object) – Two or more fitted results on the same sample, each exposing forecast_paths(steps, *, seed) returning (n_paths, steps, k); the Bayesian VAR family, the stochastic-volatility posterior, and any point-estimate result with a simulated predictive. The predictive paths must agree in steps and k across models, which is checked when paths are drawn, not here.

  • records (Sequence[_MarginalLikelihoodSelection] | None) – The models’ marginal likelihoods as MarginalLikelihoodSelection records, one per result in the same order. None computes them through marginal_likelihood(), which is exact for the conjugate BVAR and simulated (bridge sampling or Chib) for the sampled members of the family.

  • prior_probabilities (Sequence[float] | None) – \(\pi_i\), one per result, non-negative and not all zero; they are normalized to sum to one. None gives every model the same prior probability.

  • names (Sequence[str] | None) – Row labels for the record, one per result, unique. None uses each marginal-likelihood record’s source, which for a package result is the model and prior label, so two fits that differ only in prior are told apart.

Return type:

_ModelCombinationSelection

Returns:

The ModelCombinationSelection with method set to "bayesian model average", weights the posterior model probabilities, scores the log marginal likelihoods the weights came from, and n_origins zero, since no held-out evaluation was involved.

Raises:

  • SpecificationError – If fewer than two results are given; if a result exposes no forecast_paths; if a result has no marginal likelihood (a point estimate, when records is omitted); if the records score different numbers of observations; if a prior probability is negative or all are zero; or if names repeats a label.

  • DimensionError – If records, names, or prior_probabilities do not have one entry per result.

Note

When some marginal likelihoods are simulated, each carries a Monte Carlo standard error, and a weight ratio \(w_i / w_j\) is only as precise as the difference of two log marginal likelihood estimates; the record notes this when it applies. A gap of a few nats between two simulated estimates with errors of a nat each is not a decisive weight, however sharply the normalized numbers appear to split.

See also

  • stacking() – weights that maximize held-out predictive score, which do not concentrate under misspecification.

  • marginal_likelihood() – the record this function computes for each result when records is omitted.

  • compare() – the same marginal likelihoods ranked in a table, with Bayes factors, when a winner rather than a mixture is wanted.

  • ModelCombinationSelection – the record returned, and how it draws from the mixture predictive.

References

Hoeting, J. A., Madigan, D., Raftery, A. E., & Volinsky, C. T. (1999). Bayesian model averaging: A tutorial. Statistical Science, 14(4), 382-401.

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

Example

Two conjugate BVARs on one sample, differing only in the tightness of the Minnesota prior; the marginal likelihoods are exact, so the weights are too:

>>> 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 = BVAR(y, order=1).fit(n_draws=200, seed=0)
>>> tight = BVAR(y, order=1, prior=NormalInverseWishartPrior(tightness=0.02))
>>> tight = tight.fit(n_draws=200, seed=0)
>>> bma = bayesian_model_average(loose, tight, names=("loose", "tight"))
>>> bma.method, bma.weights.round(3)
('bayesian model average', array([0.998, 0.002]))
>>> bma.forecast_paths(4, seed=0).shape, bma.forecast(4).shape
((200, 4, 2), (4, 2, 3))

Prior model probabilities tilt the weights but a six-nat gap in log marginal likelihood is hard to overturn:

>>> tilted = bayesian_model_average(loose, tight, prior_probabilities=(1.0, 99.0))
>>> tilted.weights.round(3)
array([0.824, 0.176])

Models conditioning on different presamples score different samples and are refused:

>>> two = BVAR(y, order=2).fit(n_draws=100, seed=0)
>>> bayesian_model_average(loose, two)
Traceback (most recent call last):
    ...
cultivars.exceptions.SpecificationError: The marginal likelihoods score different ...