Prediction intervals for regression

What is being estimated#

For a new input $x$, a prediction interval at level $1-\alpha$ is a random set $C(x) = [\hat{l}(x), \hat{u}(x)]$ with

$$ P\big(Y \in C(X)\big) \ge 1 - \alpha $$

Where:

• $Y$ — the outcome of a fresh draw, not its conditional mean $\mathbb{E}[Y \mid X = x]$.
• The probability is over both the new pair and the training/calibration data — a marginal guarantee.

The distinction from a confidence interval for $\mathbb{E}[Y \mid X = x]$ is the irreducible term. Under a Gaussian linear model the two differ by exactly that term:

$$ \hat{y}0 \pm t{n-p,\,1-\alpha/2}\; \hat{\sigma} \sqrt{1 + x_0^\top (X^\top X)^{-1} x_0} $$

The $1$ under the root is the new observation's own noise; drop it and you have the confidence interval for the mean. As $n \to \infty$ the leverage term vanishes and the confidence interval collapses to a point, while the prediction interval converges to $\pm z_{1-\alpha/2}\hat{\sigma}$ and stops.

Quantile regression as the workhorse#

The pinball (check) loss at level $\tau$,

$$ L_\tau(y, q) = \max\big(\tau (y - q),\; (\tau - 1)(y - q)\big) $$

is minimised in expectation at the conditional $\tau$-quantile of $Y \mid X$. Fitting at $\tau = \alpha/2$ and $\tau = 1 - \alpha/2$ gives a conditionally adaptive interval with no distributional assumption — the source of the varying widths in the developer table. Koenker and Bassett (1978) established the linear case; gradient boosting and neural networks optimise the same loss directly.

Two known failure modes: independently fitted quantiles can cross, repaired by monotone rearrangement (Chernozhukov, Fernández-Val and Galichon, 2010); and the fitted quantiles are only as calibrated as the model, with no finite-sample coverage guarantee at all.

Conformalising the interval#

Split conformal on absolute residuals (Lei et al., 2018) takes $s_i = |Y_i - \hat{\mu}(X_i)|$ on a calibration set and returns $\hat{\mu}(x) \pm \hat{q}_{1-\alpha}$. Finite-sample marginal validity under exchangeability, constant width, no adaptivity — the "flat" method above.

CQR (Romano, Patterson and Candès, 2019) conformalises the quantile pair with

$$ E_i = \max\big(\hat{q}_{\alpha/2}(X_i) - Y_i,\; Y_i - \hat{q}_{1-\alpha/2}(X_i)\big) $$

and returns $[\hat{q}{\alpha/2}(x) - \hat{Q}, \; \hat{q}{1-\alpha/2}(x) + \hat{Q}]$ with $\hat{Q}$ the $\lceil (m+1)(1-\alpha) \rceil / m$ empirical quantile of ${E_i}$. This inherits validity from conformal theory and adaptivity from the quantile models — the combination is why it is the default recommendation. Locally weighted variants scale residuals by a fitted spread model (Lei and Wasserman, 2014) to reach adaptivity from a mean model instead.

Jackknife+ and CV+ (Barber, Candès, Ramdas and Tibshirani, 2021) avoid spending data on a calibration split by recycling leave-one-out or cross-validation residuals, at the cost of a $1 - 2\alpha$ worst-case guarantee (empirically much closer to $1-\alpha$).

Grading intervals properly#

Coverage and width should be combined by a proper scoring rule, not eyeballed. The interval score (Gneiting and Raftery, 2007) does exactly this:

$$ S_\alpha(l, u; y) = (u - l) + \tfrac{2}{\alpha}(l - y)\mathbb{1}{y < l} + \tfrac{2}{\alpha}(y - u)\mathbb{1}{y > u} $$

Width is paid for directly; misses are charged in proportion to how far outside they fell. Averaged pinball loss across a grid of $\tau$ gives the CRPS, the corresponding rule for a full predictive distribution, and the standard metric in probabilistic forecasting.

Distributional and Bayesian routes#

• Heteroscedastic likelihood: predict $\mu(x)$ and $\log \sigma(x)$ jointly and minimise the Gaussian negative log-likelihood. Cheap and effective; sensitive to early-training variance collapse.
• NGBoost (Duan et al., 2020) fits full parametric predictive distributions by natural-gradient boosting.
• Deep ensembles (Lakshminarayanan, Pritzel and Blundell, 2017) and MC dropout (Gal and Ghahramani, 2016) capture model uncertainty rather than noise; both are known to be miscalibrated out of the box and both are improved by a conformal wrapper.
• The split matters conceptually: aleatoric uncertainty is irreducible noise, and no amount of data narrows it; epistemic uncertainty is ignorance about the model, and data does narrow it (Kendall and Gal, 2017). A prediction interval must contain both. Ensemble spread alone contains only the second, which is the standard way ensemble-based intervals end up far too narrow.

Under shift#

Every guarantee here rests on exchangeability between calibration and deployment data. Under covariate shift with known likelihood ratios, weighted conformal restores validity (Tibshirani et al., 2019). Under drift, adaptive conformal inference (Gibbs and Candès, 2021) adjusts $\alpha$ online from realised coverage. In practice: monitor realised coverage as a live metric, sliced, and recalibrate on recent data. Recalibration needs no refit — it is one quantile of one array.

What to learn next#

• Conformal prediction — where the coverage guarantee comes from.
• Quantile regression — the model that produces the two ends.
• Proper scoring rules — how to grade coverage and width with one number.