Probability Forecasts

A probability forecast predicts how likely each outcome is instead of naming one. Formally it reports \(P(Y \in A \mid \text{information up to time } t)\) for events or ranges \(A\), and in practice it is the difference between “sales will be 500” and “20% chance under 480, 50% chance between 480 and 520, 30% chance above 520.”

The familiar examples are the weather — “70% chance of rain” against a deterministic “it will rain” — a match called at 60/25/15, or a return that has a 10% chance of falling below −2%. It is the discrete, event-flavored member of the probabilistic forecast family: a quantile forecast reports cut points and a full distribution forecast reports the whole law, while this one attaches numbers to named outcomes.

Producing them is the easy half. Logistic regression and classifiers output class probabilities, ensembles give a spread across outcomes, Bayesian methods produce posteriors directly, and a softmax or sigmoid layer emits numbers between 0 and 1. Whether those numbers mean anything is the other half.

Calibrated and useless, informative and wrong

The standard evaluation is the Brier score, \(\frac{1}{N}\sum (p_i - o_i)^2\), together with log loss and a reliability curve. The Brier score splits into three pieces — reliability, which is miscalibration and should be small; resolution, which is how much the forecast varies with the truth and should be large; and uncertainty, which belongs to the problem and cannot be changed. Splitting it apart shows why the two halves of a good forecast have to be checked separately.

import numpy as np

rng = np.random.default_rng(0)
n = 400_000
x = rng.random(n)
true_p = 0.15 + 0.7 * x ** 2
outcome = (rng.random(n) < true_p).astype(float)
base = outcome.mean()
"""Three forecasters see the same days. The truth varies with x, and the
   base rate over all days is the number a climatology forecast would give."""

sharp = np.clip(1 / (1 + np.exp(-4 * (true_p - 0.5) / 0.15)), 0.01, 0.99)
F = (("always the base rate", np.full(n, base)),
     ("the true probability", true_p),
     ("sharp but miscalibrated", sharp))

def decompose(q, y, bins=20):
    edges = np.linspace(0, 1, bins + 1)
    idx = np.clip(np.digitize(q, edges[1:-1]), 0, bins - 1)
    rel = res = 0.0
    obar = y.mean()
    for b in range(bins):
        m = idx == b
        if not m.any():
            continue
        w = m.mean()
        rel += w * (q[m].mean() - y[m].mean()) ** 2
        res += w * (y[m].mean() - obar) ** 2
    return rel, res, obar * (1 - obar)

print(f"  base rate {base:.4f} over {n:,} days")
print(f"{'forecaster':>26}{'Brier':>9}{'reliability':>13}{'resolution':>12}"
      f"{'log loss':>10}{'accuracy':>10}")
for name, q in F:
    rel, res, unc = decompose(q, outcome)
    brier = ((q - outcome) ** 2).mean()
    ll = (-np.log(np.where(outcome == 1, q, 1 - q))).mean()
    print(f"{name:>26}{brier:>9.6f}{rel:>13.6f}{res:>12.6f}{ll:>10.6f}"
          f"{np.mean((q >= 0.5) == outcome):>10.4f}")
print("  Brier = reliability - resolution + uncertainty, with uncertainty"
      f" {base * (1 - base):.6f}")


q = sharp
edges = np.linspace(0, 1, 11)
idx = np.clip(np.digitize(q, edges[1:-1]), 0, 9)
print(f"{'says about':>12}{'days':>10}{'actually happens':>19}")
for b in range(10):
    m = idx == b
    if m.sum() < 100:
        continue
    print(f"{q[m].mean():>12.4f}{m.sum():>10,}{outcome[m].mean():>19.4f}")
#   base rate 0.3822 over 400,000 days
#                 forecaster    Brier  reliability  resolution  log loss  accuracy
#       always the base rate 0.236127     0.000000    0.000000  0.665139    0.6178
#       the true probability 0.192599     0.000006    0.043363  0.569949    0.7132
#    sharp but miscalibrated 0.250721     0.054605    0.039575  1.025736    0.7132
#   Brier = reliability - resolution + uncertainty, with uncertainty 0.236127
#   says about      days   actually happens
#       0.0152   247,374             0.2381
#       0.1448    13,415             0.4309
#       0.2474     8,686             0.4633
#       0.3484     6,904             0.4820
#       0.4492     6,403             0.4926
#       0.5501     6,123             0.5092
#       0.6510     6,432             0.5247
#       0.7522     7,862             0.5413
#       0.8541    11,396             0.5693
#       0.9790    85,405             0.7079

What the table settles

Perfect calibration is not the goal, because the forecast that always reports the base rate achieves it and says nothing. Its reliability is zero and its resolution is zero, and its Brier score is exactly the uncertainty of the problem. Any forecast worth making has to beat that, and beating it means having resolution — saying different things on different days, correctly.

The third row is the case worth remembering. It has almost as much resolution as the truth, so it genuinely knows which days are riskier, and its numbers are still not probabilities: days it calls 0.0152 happen 23.81% of the time. Its miscalibration outweighs its information, and it scores worse than the forecast that says the same thing every day. Accuracy ranks it level with the truth, because both cross 0.5 on the same days. A model that is useful for ranking can be unusable for any decision that reads the number as a probability, and only the decomposition shows it — which is what recalibration exists to repair.

What they are good for

Probabilities are usable directly by a decision rule in a way that a class label is not. A label answers “rain or not”; a probability combines with the cost of carrying an umbrella against the cost of being soaked, and the threshold falls out of that ratio rather than sitting at 0.5. That is the whole advantage, and it survives only if the numbers are calibrated.

They are also the most interpretable output a model produces. “70% chance of rain” needs no explanation, which is why the format spread from weather to finance, sport and risk management. The interpretability is a liability when the calibration is poor, since nobody reading the number knows to discount it.

Inspect calibration with reliability curves and compare proper scores against the base-rate forecast on held-out data. A worse uncalibrated score can coexist with useful ranking or resolution, as the example shows. Fit a recalibration map on separate data and assess its held-out score before discarding the model; deployment against probability thresholds requires validation at the operating range.

References

• NumPy documentation, numpy.random.Generator and numpy.digitize. numpy.org