3. Conformal Prediction with a Point Forecaster

Why split conformal fits a point forecaster

A standard regression model — gradient boosting, random forest, ridge regression, or a small neural network — outputs a single number ŷ(x). It gives no honest statement about how far that number might be from the truth. Split conformal prediction turns any such point forecaster into an interval generator with a finite-sample, distribution-free coverage guarantee, provided only that the calibration and test points are exchangeable.

The recipe is deliberately simple. Fit the model on a proper training split. Compute absolute residuals on a disjoint calibration split. Take the appropriate empirical quantile of those residuals and add/subtract it from the point prediction at test time. The only assumption is exchangeability; no normality, homoscedasticity, or model correctness is required.

Formal construction

Let the training set be split into D_train (size n) and D_cal (size m). Fit f̂ on D_train. For each calibration point (x_i, y_i), define the nonconformity score

s_i = |y_i − f̂(x_i)|,   i = 1, …, m

For a target miscoverage level α (e.g. 0.10 for 90% coverage), compute

q̂ = Quantile( {s_1, …, s_m}, ⌈(m+1)(1−α)⌉ / m )

using the higher empirical quantile (the ⌈·⌉ index, not interpolation). The prediction interval at a new point x is

C(x) = [ f̂(x) − q̂ ,  f̂(x) + q̂ ]

Under exchangeability of the m+1 calibration-plus-test points, P(Y ∈ C(X)) ≥ 1 − α exactly, for any f̂. The guarantee is marginal (averaged over calibration and test draws), not conditional on x.

Reference implementation

The following is a self-contained, runnable script. It uses a synthetic heteroscedastic regression problem so that interval width should visibly track the noise level, and it reports empirical coverage and mean width on both calibration and a fresh test split.

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(0)

def make_data(n):
    X = rng.uniform(-3, 3, size=(n, 2))
    noise = 0.3 + 0.5 * np.abs(X[:, 0])          # heteroscedastic
    y = np.sin(X[:, 0]) + 0.5 * X[:, 1] + rng.normal(0, noise)
    return X, y

X, y = make_data(4000)
X_tr, X_tmp, y_tr, y_tmp = train_test_split(X, y, test_size=0.5, random_state=0)
X_cal, X_te, y_cal, y_te = train_test_split(X_tmp, y_tmp, test_size=0.5, random_state=0)

model = GradientBoostingRegressor(n_estimators=300, max_depth=3, random_state=0)
model.fit(X_tr, y_tr)

alpha = 0.10
m = len(y_cal)
scores = np.abs(y_cal - model.predict(X_cal))
k = int(np.ceil((m + 1) * (1 - alpha)))
q_hat = np.sort(scores)[min(k, m) - 1]           # higher empirical quantile

def evaluate(Xs, ys, name):
    pred = model.predict(Xs)
    lo, hi = pred - q_hat, pred + q_hat
    cov = np.mean((ys >= lo) & (ys <= hi))
    width = np.mean(hi - lo)
    print(f"{name:12s} n={len(ys):4d}  coverage={cov:.3f}  mean_width={width:.3f}")

print(f"q_hat = {q_hat:.4f}  (alpha={alpha}, m={m})")
evaluate(X_cal, y_cal, "calibration")
evaluate(X_te,  y_te,  "test")

Expected behaviour on this synthetic setup: calibration coverage sits at or just above 1 − α (the finite-sample guarantee is conservative on the calibration set itself), and test coverage lands within roughly ±0.02 of 0.90 for m ≈ 1000. Mean width is a single scalar 2·q̂ — the interval is constant-width in ŷ, which is the main efficiency limitation of split conformal with absolute residuals.

Reading the calibration–test relationship

The diagram below shows how the calibration residuals collapse into a single threshold q̂, and how that threshold is then applied uniformly to test predictions. It also highlights the two failure modes worth checking: a calibration set too small to resolve the target quantile, and a test distribution that has drifted away from the calibration distribution.

Diagnostics that matter in practice

Two numbers decide whether the wrapper is doing its job: empirical coverage on the test split and mean interval width. Coverage below 1 − α by more than sampling noise (roughly 1.96·√(α(1−α)/n_test)) signals a broken assumption — most often leakage between D_train and D_cal, or a test set drawn from a shifted distribution. Width that is far larger than the residual spread suggests m is too small: with m = 100 and α = 0.10, the quantile index is ⌈101·0.9⌉ = 91, so q̂ is the 91st of 100 residuals and is highly variable across calibration draws.

A useful sanity check is to repeat the calibration split several times with different seeds and record the spread of q̂. If the standard deviation of q̂ exceeds a few percent of its mean, increase m or move to a cross-conformal scheme that reuses all training data.

What to hand off

The deliverable from this case is a reusable wrapper: a function that takes (X_train, y_train, X_cal, y_cal, X_test, alpha), fits the point model, computes q̂ with the correct higher quantile, and returns (lower, upper, q_hat). Pair it with the coverage/width report above so that every downstream use of the intervals is accompanied by the empirical evidence that the guarantee is holding on the data at hand.