Skip to content

Stacking

Bayesian model stacking and ensemble prediction.

trade_study.stack_bayesian(compare_dict, *, method='stacking')

Bayesian stacking via arviz.compare.

Parameters:

Name Type Description Default
compare_dict dict[str, Any]

Dictionary mapping model names to arviz DataTree or ELPDData objects (must contain log_likelihood group).

 required
method str

Weighting method. One of "stacking", "BB-pseudo-BMA", "pseudo-BMA".

 'stacking'

Returns:

Type Description
dict[str, float]

Dictionary mapping model names to stacking weights.
Source code in src/trade_study/stacking.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def stack_bayesian(
    compare_dict: dict[str, Any],
    *,
    method: str = "stacking",
) -> dict[str, float]:
    """Bayesian stacking via arviz.compare.

    Args:
        compare_dict: Dictionary mapping model names to arviz DataTree
            or ELPDData objects (must contain log_likelihood group).
        method: Weighting method. One of "stacking", "BB-pseudo-BMA",
            "pseudo-BMA".

    Returns:
        Dictionary mapping model names to stacking weights.
    """
    import arviz as az  # type: ignore[import-untyped]

    result = az.compare(compare_dict, method=method)
    return dict(zip(result.index, result["weight"], strict=True))

trade_study.stack_scores(score_matrix, *, maximize=False)

Optimize stacking weights from a score matrix.

For non-Bayesian models where log-likelihoods aren't available. The objective is linear in the supplied model scores, so the solution puts weight only on models with the best mean score. Exactly equal computed best means split weight uniformly; near ties are not treated as equal. Positive changes of score units preserve the optimum, subject to floating precision. No iterative optimizer or SciPy dependency is needed. This does not optimize the score of a mixture distribution or the loss of an averaged prediction.

Parameters:

Name Type Description Default
score_matrix NDArray[floating[Any]]

Non-empty finite array of shape (n_models, n_test_points). Each entry is the score of model i on test point j.

 required
maximize bool

If True, maximize the weighted score; if False, minimize.

 False

Returns:

Type Description
NDArray[floating[Any]]

Non-negative float64 weights, shape (n_models,), summing to 1.

Raises:

Type Description
ValueError

If scores are empty, non-finite, or not two-dimensional.
FloatingPointError

If model means cannot be represented numerically.

Examples:

>>> weights = stack_scores(np.array([[1.0, 1.0], [3.0, 3.0]]))
>>> weights.tolist()
[1.0, 0.0]
>>> stack_scores(1e6 * np.array([[1.0, 1.0], [3.0, 3.0]])).tolist()
[1.0, 0.0]
>>> stack_scores(np.array([[1.0], [1.0], [3.0]])).tolist()
[0.5, 0.5, 0.0]
Source code in src/trade_study/stacking.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def stack_scores(
    score_matrix: NDArray[np.floating[Any]],
    *,
    maximize: bool = False,
) -> NDArray[np.floating[Any]]:
    """Optimize stacking weights from a score matrix.

    For non-Bayesian models where log-likelihoods aren't available.
    The objective is linear in the supplied model scores, so the solution puts
    weight only on models with the best mean score. Exactly equal computed
    best means split weight uniformly; near ties are not treated as equal.
    Positive changes of score units preserve the optimum, subject to floating
    precision. No iterative optimizer or SciPy dependency is needed.
    This does not optimize the score of a mixture distribution or the loss
    of an averaged prediction.

    Args:
        score_matrix: Non-empty finite array of shape (n_models, n_test_points). Each
            entry is the score of model i on test point j.
        maximize: If True, maximize the weighted score; if False, minimize.

    Returns:
        Non-negative float64 weights, shape (n_models,), summing to 1.

    Raises:
        ValueError: If scores are empty, non-finite, or not two-dimensional.
        FloatingPointError: If model means cannot be represented numerically.

    Examples:
        >>> weights = stack_scores(np.array([[1.0, 1.0], [3.0, 3.0]]))
        >>> weights.tolist()
        [1.0, 0.0]
        >>> stack_scores(1e6 * np.array([[1.0, 1.0], [3.0, 3.0]])).tolist()
        [1.0, 0.0]
        >>> stack_scores(np.array([[1.0], [1.0], [3.0]])).tolist()
        [0.5, 0.5, 0.0]
    """
    scores = np.asarray(score_matrix, dtype=np.float64)
    if scores.ndim != 2 or 0 in scores.shape or not np.all(np.isfinite(scores)):
        msg = "score_matrix must be a non-empty, finite two-dimensional array"
        raise ValueError(msg)
    # Divide before summing and use the extended accumulator where available.
    # An unrepresentable average must raise instead of producing unchecked weights.
    try:
        with np.errstate(over="raise", invalid="raise"):
            means = np.sum(scores.astype(np.longdouble) / scores.shape[1], axis=1)
    except FloatingPointError as error:
        msg = "score_matrix means are not representable; rescale the score units"
        raise FloatingPointError(msg) from error
    best = np.max(means) if maximize else np.min(means)
    selected = means == best
    return np.asarray(selected / np.count_nonzero(selected), dtype=np.float64)

trade_study.stack_proportional(score_matrix, *, maximize=False, temperature=1.0)

Weight models proportionally to exponential score utilities.

Mean scores are transformed using a direction-aware softmax. Near-tied models receive comparable weights in either direction, including when scores are zero or negative. These are heuristic ensemble weights, rather than Bayesian stacking weights or calibrated probabilities.

Parameters:

Name Type Description Default
score_matrix NDArray[floating[Any]]

Finite array of shape (n_models, n_test_points).

 required
maximize bool

Whether higher scores are better. Default is minimize.

 False
temperature float

Positive, finite scale in score units. Larger values spread weight more evenly; smaller values favor the best model. The default is 1.0. Scale this with the scores when changing their units to preserve the same weights.

 1.0

Returns:

Type Description
NDArray[floating[Any]]

Non-negative float64 weights summing to 1, shape (n_models,).
NDArray[floating[Any]]

Exactly equal mean scores receive uniform weights.

Raises:

Type Description
ValueError

If scores are empty, non-finite, or not two-dimensional, or temperature is not positive and finite.
Source code in src/trade_study/stacking.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def stack_proportional(
    score_matrix: NDArray[np.floating[Any]],
    *,
    maximize: bool = False,
    temperature: float = 1.0,
) -> NDArray[np.floating[Any]]:
    """Weight models proportionally to exponential score utilities.

    Mean scores are transformed using a direction-aware softmax. Near-tied
    models receive comparable weights in either direction, including when
    scores are zero or negative. These are heuristic ensemble weights,
    rather than Bayesian stacking weights or calibrated probabilities.

    Args:
        score_matrix: Finite array of shape (n_models, n_test_points).
        maximize: Whether higher scores are better. Default is minimize.
        temperature: Positive, finite scale in score units. Larger values
            spread weight more evenly; smaller values favor the best model.
            The default is 1.0. Scale this with the scores when changing
            their units to preserve the same weights.

    Returns:
        Non-negative float64 weights summing to 1, shape (n_models,).
        Exactly equal mean scores receive uniform weights.

    Raises:
        ValueError: If scores are empty, non-finite, or not two-dimensional,
            or temperature is not positive and finite.
    """
    scores = np.asarray(score_matrix, dtype=np.float64)
    if scores.ndim != 2 or 0 in scores.shape or not np.all(np.isfinite(scores)):
        msg = "score_matrix must be a non-empty, finite two-dimensional array"
        raise ValueError(msg)
    if not np.isfinite(temperature) or temperature <= 0:
        msg = "temperature must be positive and finite"
        raise ValueError(msg)

    # Divide before summing to avoid overflowing the mean of large scores.
    means = np.sum(scores / scores.shape[1], axis=1)
    utility = means if maximize else -means
    # Overflow in a very large unfavorable gap means a zero softmax weight.
    with np.errstate(over="ignore"):
        logits = (utility - np.max(utility)) / temperature
    weights = np.exp(logits)
    return np.asarray(weights / weights.sum(), dtype=np.float64)

trade_study.ensemble_predict(predictions, weights)

Weighted ensemble of model predictions.

Parameters:

Name Type Description Default
predictions list[NDArray[floating[Any]]]

List of prediction arrays, one per model. Each should have the same shape.

 required
weights NDArray[floating[Any]]

Stacking weights, shape (n_models,). Normalization leaves the supplied array and its backing storage unchanged. Read-only inputs and views are supported.

 required

Returns:

Type Description
NDArray[floating[Any]]

Weighted average prediction array.

Examples:

>>> import numpy as np
>>> predictions = [np.array([1.0, 2.0]), np.array([3.0, 4.0])]
>>> ensemble_predict(predictions, np.array([0.25, 0.75])).tolist()
[2.5, 3.5]
>>> weights = np.array([1.0, 3.0])
>>> weights.setflags(write=False)
>>> ensemble_predict(predictions, weights).tolist(), weights.tolist()
([2.5, 3.5], [1.0, 3.0])
Source code in src/trade_study/stacking.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def ensemble_predict(
    predictions: list[NDArray[np.floating[Any]]],
    weights: NDArray[np.floating[Any]],
) -> NDArray[np.floating[Any]]:
    """Weighted ensemble of model predictions.

    Args:
        predictions: List of prediction arrays, one per model.
            Each should have the same shape.
        weights: Stacking weights, shape (n_models,). Normalization leaves the
            supplied array and its backing storage unchanged. Read-only inputs
            and views are supported.

    Returns:
        Weighted average prediction array.

    Examples:
        >>> import numpy as np
        >>> predictions = [np.array([1.0, 2.0]), np.array([3.0, 4.0])]
        >>> ensemble_predict(predictions, np.array([0.25, 0.75])).tolist()
        [2.5, 3.5]
        >>> weights = np.array([1.0, 3.0])
        >>> weights.setflags(write=False)
        >>> ensemble_predict(predictions, weights).tolist(), weights.tolist()
        ([2.5, 3.5], [1.0, 3.0])
    """
    w = np.array(weights, dtype=np.float64, copy=True)
    w /= w.sum()
    result = np.zeros_like(predictions[0], dtype=np.float64)
    for pred, wi in zip(predictions, w, strict=True):
        result += wi * pred
    return result