4. Jackknife+ Prediction Intervals

Why Jackknife+ Exists

Split conformal prediction is cheap and distribution-free, but it pays a price: the model is trained on only part of the data, and the calibration set is disjoint from the training set. With small n, that split wastes information and produces intervals that are wider than necessary. The jackknife+ construction removes the split entirely. It trains n leave-one-out (LOO) models, uses each held-out point as its own calibration residual, and aggregates the resulting intervals in a way that preserves a rigorous finite-sample coverage guarantee — specifically, at least 1 − 2α marginal coverage for a nominal 1 − α target.

The key subtlety is that jackknife+ does not achieve 1 − α coverage; it achieves 1 − 2α. This is the price of using LOO residuals, which are not exchangeable with the test residual in the same way split-conformal calibration residuals are. The jackknife+ after bootstrap (JaB) variant replaces the n LOO fits with B bootstrap fits, trading a small additional coverage slack for far lower computational cost when n is large.

Formal Construction

Let the training data be D = {(x_i, y_i)}_{i=1..n} and let μ̂_{−i} denote the model trained on D \ {(x_i, y_i)}. Define the LOO residuals R_i = |y_i − μ̂_{−i}(x_i)|. For a test point x, jackknife+ forms the interval

[ q⁻_{α}( μ̂_{−i}(x) − R_i ),  q⁺_{α}( μ̂_{−i}(x) + R_i ) ]

where q⁻ and q⁺ are the α and 1−α empirical quantiles over the n values {μ̂_{−i}(x) − R_i} and {μ̂_{−i}(x) + R_i} respectively. The lower endpoint is the α-quantile of the lower candidate set, and the upper endpoint is the (1−α)-quantile of the upper candidate set. This asymmetric quantile pairing is what distinguishes jackknife+ from a naive "residual quantile" interval and is what makes the coverage guarantee hold.

For JaB, draw B bootstrap resamples D^{*b}, fit μ̂^{*b}, and compute out-of-bag residuals R_i^{*b} for each observation i not in D^{*b}. The interval endpoints become quantiles over the pooled set of {μ̂^{*b}(x) ± R_i^{*b}} pairs. The guarantee weakens to 1 − 2α − O(1/√B) in the worst case, but in practice the slack is small for B ≥ 200.

Reference Implementation

The following Python implementation uses scikit-learn's LinearRegression as the base model. It is self-contained: it generates synthetic data, fits jackknife+, jackknife+ after bootstrap, and split conformal, then reports empirical coverage and mean interval width on a held-out test set.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(42)

def make_data(n, p=5, noise=1.0):
    X = rng.normal(size=(n, p))
    beta = rng.normal(size=p)
    y = X @ beta + noise * rng.normal(size=n)
    return X, y

def fit_predict(X_tr, y_tr, X_te):
    m = LinearRegression().fit(X_tr, y_tr)
    return m.predict(X_te)

def jackknife_plus(X_tr, y_tr, X_te, alpha=0.1):
    n = len(y_tr)
    lo = np.empty((n, len(X_te)))
    hi = np.empty((n, len(X_te)))
    for i in range(n):
        mask = np.ones(n, dtype=bool); mask[i] = False
        pred_te = fit_predict(X_tr[mask], y_tr[mask], X_te)
        pred_i  = fit_predict(X_tr[mask], y_tr[mask], X_tr[i:i+1])[0]
        R_i = abs(y_tr[i] - pred_i)
        lo[i] = pred_te - R_i
        hi[i] = pred_te + R_i
    lower = np.quantile(lo, alpha, axis=0)
    upper = np.quantile(hi, 1 - alpha, axis=0)
    return lower, upper

def jab(X_tr, y_tr, X_te, alpha=0.1, B=200):
    n = len(y_tr)
    lo_all, hi_all = [], []
    for b in range(B):
        idx = rng.integers(0, n, size=n)
        oob = np.setdiff1d(np.arange(n), idx)
        if len(oob) == 0:
            continue
        pred_te = fit_predict(X_tr[idx], y_tr[idx], X_te)
        pred_oob = fit_predict(X_tr[idx], y_tr[idx], X_tr[oob])
        R = np.abs(y_tr[oob] - pred_oob)
        lo_all.append(pred_te[None, :] - R[:, None])
        hi_all.append(pred_te[None, :] + R[:, None])
    lo = np.vstack(lo_all); hi = np.vstack(hi_all)
    return np.quantile(lo, alpha, axis=0), np.quantile(hi, 1 - alpha, axis=0)

def split_conformal(X_tr, y_tr, X_te, alpha=0.1):
    X_fit, X_cal, y_fit, y_cal = train_test_split(
        X_tr, y_tr, test_size=0.3, random_state=0)
    pred_cal = fit_predict(X_fit, y_fit, X_cal)
    R = np.abs(y_cal - pred_cal)
    q = np.quantile(R, 1 - alpha)
    pred_te = fit_predict(X_fit, y_fit, X_te)
    return pred_te - q, pred_te + q

def evaluate(lo, hi, y_te):
    cov = np.mean((y_te >= lo) & (y_te <= hi))
    width = np.mean(hi - lo)
    return cov, width

X, y = make_data(n=300)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=500, random_state=1)

for name, fn in [("Jackknife+", jackknife_plus),
                 ("JaB (B=200)", jab),
                 ("Split conformal", split_conformal)]:
    lo, hi = fn(X_tr, y_tr, X_te, alpha=0.1)
    cov, w = evaluate(lo, hi, y_te)
    print(f"{name:18s} coverage={cov:.3f}  width={w:.3f}")

Expected output shape (values will vary with the random seed, but the ordering is stable): jackknife+ and JaB should show coverage near or above 0.90 on the test set, with widths slightly narrower than split conformal because they use all n training points. Split conformal will typically land near 0.90 as well but with wider intervals, especially when n is small. If you observe jackknife+ coverage dropping well below 0.80, the most common causes are (a) a non-exchangeable data split (e.g., time-ordered data), (b) a bug where the LOO residual uses the full-data prediction instead of the held-out prediction, or (c) an incorrect quantile level — remember the lower endpoint uses α, not 1−α.

Comparing the Three Methods

The practical decision rule: if n is small (say under a few thousand) and the base model is cheap to refit, use jackknife+ for the tightest intervals. If n is large or the model is expensive (e.g., gradient boosting with many trees), use JaB with B = 200–500; the coverage slack is negligible in practice and the runtime drops by orders of magnitude. Split conformal remains the right choice when you need the exact 1 − α guarantee and can afford to sacrifice training data.