stacking#

cultivars.bayes.combination.stacking(log_density, results, *, names=None)[source]#

Weight models to maximize the mixture’s held-out log predictive score.

Given the log predictive density \(\ell_{tm} = \log p_m(y_t \mid y_{<t})\) that each model \(m\) assigned to the realized outcome at each evaluation origin \(t\), the weights solve

\[\max_{w \in \Delta^{M-1}} \; \frac{1}{T} \sum_{t=1}^{T} \log \sum_{m=1}^{M} w_m\, e^{\ell_{tm}},\]

the mean log score of the mixture of predictive densities over the simplex (Yao, Vehtari, Simpson & Gelman, 2018). The objective is concave in \(w\), and it is maximized through a softmax map from \(M - 1\) free coordinates whose stationary point is the simplex optimum, so the weights are non-negative and sum to one by construction rather than by projection.

Nothing here assumes a model is true. Two models that err in different directions keep sharing weight however long the evaluation window, which is why stacking is the combination for forecasting and why a weight vector that does not concentrate is a finding about the models, not a failure to choose. The price is that the densities must be genuinely held out: from a rolling or expanding origin, each \(\ell_{tm}\) computed from a fit that did not see \(y_t\). In-sample densities reward whichever model overfits most.

Parameters:

  • log_density (ArrayLike) – (T, M) log predictive densities of the realized outcome, one row per evaluation origin and one column per model, in the order of results; anything numpy.asarray() accepts, finite throughout. A rolling backtest produces exactly this matrix. At least M origins are required, though a useful weight needs many more; the record reports the count as n_origins.

  • results (Sequence[object]) – The M fitted results the columns correspond to, each exposing forecast_paths(steps, *, seed); typically the fits on the full sample, which the mixture predicts from.

  • names (Sequence[str] | None) – Row labels for the record, one per result, unique. None uses model[0] ... model[M-1].

Return type:

_ModelCombinationSelection

Returns:

The ModelCombinationSelection with method set to "stacking", weights the optimal simplex point, scores each model’s own mean held-out log density (column means of log_density), and n_origins the number of rows.

Raises:

  • DimensionError – If log_density is not two-dimensional, if its column count differs from the number of results, or if names does not have one label per result.

  • SpecificationError – If there are fewer than two models, fewer origins than models, a result exposes no forecast_paths, or names repeats a label.

  • NumericalError – If a log density is not finite, or if the optimizer fails to converge to a finite objective.

Note

The record’s first note reports the stacked mean log score against the best single model’s. The stacked score can never be below the best single model’s on the evaluation origins, since a vertex of the simplex is feasible; how much it exceeds it is the gain from combining, on the data the weights were chosen on, and is therefore optimistic.

See also

  • bayesian_model_average() – weights from marginal likelihoods, for the case where one model is believed true.

  • Backtest – the rolling-origin machinery that produces held-out predictive densities.

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

References

Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using stacking to average Bayesian predictive distributions. Bayesian Analysis, 13(3), 917-1007.

Geweke, J., & Amisano, G. (2011). Optimal prediction pools. Journal of Econometrics, 164(1), 130-141.

Example

Two stand-in models – any object with forecast_paths will do – and held-out log densities in which the first is better at every origin, so the optimum is the vertex:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> class Toy:
...     def __init__(self, level):
...         self.level = level
...     def forecast_paths(self, steps=8, *, seed=None):
...         return np.full((100, steps, 1), self.level)
>>> good = -0.5 * rng.standard_normal(200) ** 2
>>> logdens = np.column_stack([good, good - 0.3])
>>> combined = stacking(logdens, (Toy(0.0), Toy(1.0)), names=("a", "b"))
>>> combined.weights.round(3), combined.n_origins
(array([1., 0.]), 200)
>>> float(combined.forecast_paths(4, seed=0).mean())
0.0

Two models whose errors are unrelated share the weight, and that is the point:

>>> other = -0.5 * rng.standard_normal(200) ** 2
>>> shared = stacking(np.column_stack([good, other]), (Toy(0.0), Toy(1.0)))
>>> shared.names, bool(np.all((shared.weights > 0.3) & (shared.weights < 0.7)))
(('model[0]', 'model[1]'), True)
On this page