Prediction Intervals and Split Conformal Prediction › Implementation in NumPy and scikit-learn
Machine Learning from ScratchAll chapters Prediction Intervals and Split Conformal Prediction › Implementation in NumPy and scikit-learnSection view · Sections in this chapter6 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
We build the conformal quantile and interval arithmetic from NumPy primitives, using scikit-learn for the base predictors. Fitting and calibration data construct the intervals; untouched test data evaluate them. The example keeps all three roles disjoint and performs no model or score selection using calibration or test outcomes.
The data are one-dimensional and heteroscedastic, with the noise standard deviation growing quadratically away from the origin. This is a realistic situation where a constant-width interval will struggle, so it doubles as a stress test.
PythonIn [7]Expand
def make_data(n, rng, shift=0.0):
"""y = sin(x) + 0.5x + noise whose sd grows with |x|."""
x = rng.uniform(-3.0, 3.0, size=n) + shift
sigma = 0.30 + 0.25 * x**2 # heteroscedastic noise scale
y = np.sin(x) + 0.5 * x + rng.normal(0.0, 1.0, size=n) * sigma
return x, y
x_all, y_all = make_data(6000, rng)
# Three disjoint roles: fit, calibrate, evaluate.
x_tr, y_tr = x_all[:3000], y_all[:3000]
x_ca, y_ca = x_all[3000:4500], y_all[3000:4500]
x_te, y_te = x_all[4500:], y_all[4500:]We fit an ordinary least-squares line even though the synthetic mean is nonlinear. This makes model misspecification visible. The held-out calibration procedure still has its marginal coverage guarantee under the sampling assumptions, but the predictor and score can produce inefficient intervals or uneven local coverage.
PythonIn [8]Expand
X_tr, X_ca, X_te = (a.reshape(-1, 1) for a in (x_tr, x_ca, x_te))
model = LinearRegression().fit(X_tr, y_tr)
pred_ca = model.predict(X_ca)
pred_te = model.predict(X_te)conformal_quantile implements the explicit order-statistic rule and its infinite-threshold case. regression_intervals adds the calibrated margin to predictions. Sorting and indexing make the finite-sample convention visible; calling an ordinary interpolating quantile at is not generally equivalent.
PythonIn [9]Expand
def conformal_quantile(scores, alpha):
"""Finite-sample conformal quantile: the ceil((n+1)(1-alpha))-th smallest score."""
scores = np.asarray(scores, dtype=float)
if scores.ndim != 1 or scores.size == 0 or not np.all(np.isfinite(scores)):
raise ValueError(
"scores must be a nonempty finite one-dimensional array"
)
if not np.isscalar(alpha) or not np.isfinite(alpha) or not 0 < alpha < 1:
raise ValueError("alpha must lie strictly between zero and one")
n = scores.size
k = int(np.ceil((n + 1) * (1 - alpha)))
if k > n:
return np.inf, k
return float(np.sort(scores)[k - 1]), k
def regression_intervals(pred, q_hat, scale=1.0):
"""Symmetric conformal band; `scale` lets us plug in a local noise estimate."""
return pred - q_hat * scale, pred + q_hat * scale
# Small calibration samples must use the full-space branch rather than clipping.
assert np.isinf(conformal_quantile([0.5, 0.8, 1.2, 2.5], 0.10)[0])
assert conformal_quantile([0.5, 0.8, 1.2, 2.5], 0.20) == (2.5, 4)
assert conformal_quantile(np.arange(19), 0.10) == (17.0, 18)With those in hand, the entire method is six lines.
PythonIn [10]Expand
alpha = 0.10
scores_ca = np.abs(y_ca - pred_ca)
q_hat, k_ca = conformal_quantile(scores_ca, alpha)
lo, hi = regression_intervals(pred_te, q_hat)
covered = (y_te >= lo) & (y_te <= hi)ConsoleOut [11]Expand
n_cal = 1500, k = 1351, q_hat = 2.338 empirical test coverage = 0.909 (nominal target = 0.90) mean interval width = 4.675
The printed coverage describes this realized held-out sample. Constant-width intervals can cover more often near the low-noise center and less often in the noisy tails. The mean width alone does not reveal those local differences, and a result near nominal in this run is not proof of conditional or per-sample validity.
The sweep applies several nominal levels to the same calibration scores and test data. Each row changes the quantile cutoff and recomputes inclusion; for levels that require the infinite-threshold branch the row reports the full response space and its empirical coverage is one. The rows are related evaluations of one sample, not independent replications or a guarantee that every measured coverage exceeds nominal.
PythonIn [12]Expand
alphas = np.array([0.02, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50])
rows = []
for a in alphas:
q, kk = conformal_quantile(scores_ca, a)
l, h = regression_intervals(pred_te, q)
rows.append(
{
"alpha": a,
"nominal": 1 - a,
"k": kk,
"q_hat": q,
"test_coverage": ((y_te >= l) & (y_te <= h)).mean(),
"mean_width": np.mean(h - l),
}
)
coverage_table = pd.DataFrame(rows)| alpha | nominal | k | q_hat | test_coverage | mean_width | |
|---|---|---|---|---|---|---|
| 0 | 0.02 | 0.98 | 1471 | 3.859 | 0.974 | 7.719 |
| 1 | 0.05 | 0.95 | 1426 | 3.083 | 0.957 | 6.167 |
| 2 | 0.10 | 0.90 | 1351 | 2.338 | 0.909 | 4.675 |
| 3 | 0.20 | 0.80 | 1201 | 1.550 | 0.817 | 3.100 |
| 4 | 0.30 | 0.70 | 1051 | 1.086 | 0.720 | 2.172 |
| 5 | 0.40 | 0.60 | 901 | 0.818 | 0.619 | 1.636 |
| 6 | 0.50 | 0.50 | 751 | 0.634 | 0.524 | 1.267 |
The interval width grows as nominal coverage rises. Empirical test coverage fluctuates around the target and can fall below it for a particular calibration sample and finite test set; the marginal theorem does not require every row of this table to exceed nominal.
Constants and Calibration Size
Three key design choices are the miscoverage level , the calibration size , and the score. Understanding their roles is what separates a working deployment from a mysterious one.
The level specifies the tolerated marginal miscoverage. Its choice depends on the decision, the cost of misses and the usefulness of the resulting sets. Calibration size constrains whether that level permits a finite cutoff, and the experiment uses as an illustrative choice.
The calibration size controls rank granularity. With distinct scores the coverage excess from rounding is less than ; this is a probability-scale bound, not a universal bound on excess interval width. For example, and give . Quantile variability also depends on the score distribution, especially its density near the cutoff, so there is no distribution-free sample count that guarantees a stable narrow interval.
The score is the main modelling lever. Absolute residuals give constant width. Normalized residuals give locally adaptive width. Signed residuals give one-sided bounds. We turn to the adaptive case next, because it is where the constant-width weakness becomes visible.
Normalized Scores
Absolute-residual scores use the same cutoff at every input. An input-dependent positive scale can redistribute width. Fix that scale using fitting data only, then calibrate . The same exchangeability proof applies, because the score is fixed before calibration. Here is an error-scale proxy, not a validated conditional standard deviation; the construction does not guarantee improved conditional coverage.
We build by regressing the training absolute residuals on , then clip it away from zero to avoid division blow-ups. This is a crude but honest noise model, and it is fitted only on training data, consistent with the requirement that the score function be frozen before calibration.
PythonIn [14]Expand
abs_res_tr = np.abs(y_tr - model.predict(X_tr))
scale_model = LinearRegression().fit(np.abs(x_tr).reshape(-1, 1), abs_res_tr)
def sigma_hat(x):
"""Local noise scale, fitted on TRAIN residuals only."""
return np.maximum(scale_model.predict(np.abs(x).reshape(-1, 1)), 1e-2)
scores_ca_n = np.abs(y_ca - pred_ca) / sigma_hat(x_ca)
q_norm, k_norm = conformal_quantile(scores_ca_n, alpha)
lo_n, hi_n = regression_intervals(pred_te, q_norm, sigma_hat(x_te))
covered_n = (y_te >= lo_n) & (y_te <= hi_n)ConsoleOut [15]Expand
constant-width coverage = 0.909, mean width = 4.675 normalized-score coverage = 0.898, mean width = 4.138
The printed results show how normalization changes both width and measured coverage. The normalized version allocates less width near the origin and more width at larger , but this crude scale model can still under-cover a region. Its empirical test coverage can also lie below nominal. Both score choices have the same type of marginal guarantee, rather than identical realized coverage or guaranteed local accuracy.
Next we demonstrate the failure mode. We keep the model and the calibrated quantile fixed, then evaluate on data drawn from a shifted input distribution. The input shift breaks exchangeability with calibration, so the ordinary coverage theorem no longer applies. The following evaluation measures what happens in this sample.
PythonIn [16]Expand
x_sh, y_sh = make_data(
1500, np.random.default_rng(99), shift=1.5
) # inputs pushed right
pred_sh = model.predict(x_sh.reshape(-1, 1))
lo_sh, hi_sh = regression_intervals(pred_sh, q_hat)
cov_sh = ((y_sh >= lo_sh) & (y_sh <= hi_sh)).mean()ConsoleOut [17]Expand
exchangeable test coverage = 0.909 shifted test coverage = 0.807 (nominal target = 0.90)
The shifted set is covered far less often than nominal. Nothing about the model or the quantile changed; only the assumption did. The guarantee depends on the specified set-valued procedure and the joint distribution of calibration and test scores. When the deployment distribution differs from the calibration distribution, coverage can fall below nominal and must be monitored empirically.
Classification follows the same recipe with a different score. We generate a three-class problem, fit a multinomial logistic regression, and conformalize it.
PythonIn [18]Expand
Xc, yc = make_classification(
n_samples=6000,
n_features=8,
n_informative=5,
n_redundant=0,
n_classes=3,
n_clusters_per_class=1,
class_sep=0.9,
random_state=0,
)
Xc_tr, yc_tr = Xc[:3000], yc[:3000]
Xc_ca, yc_ca = Xc[3000:4500], yc[3000:4500]
Xc_te, yc_te = Xc[4500:], yc[4500:]
clf = LogisticRegression(max_iter=1000).fit(Xc_tr, yc_tr)
p_ca = clf.predict_proba(Xc_ca)
p_te = clf.predict_proba(Xc_te)
# Nonconformity score: 1 minus the probability assigned to the true label.
s_cls = 1.0 - p_ca[np.arange(yc_ca.size), yc_ca]
q_cls, k_cls = conformal_quantile(s_cls, alpha)
threshold = 1.0 - q_cls
in_set = p_te >= threshold # (n_test, n_classes) boolean
set_sizes = in_set.sum(axis=1)
covered_cls = in_set[np.arange(yc_te.size), yc_te]ConsoleOut [19]Expand
n_cal = 1500, k = 1351, q_hat = 0.690 probability threshold = 0.310 empirical test coverage = 0.895 (nominal target = 0.90) mean prediction set size = 1.18 singleton sets = 0.819 empty sets = 0.000
The printed classification coverage can be slightly below nominal on this finite test sample. In this seeded run there are no empty sets: most inputs receive one label, some receive two, and one receives all three. A calibration-set check is different from a test-set guarantee: by construction, the fraction of calibration scores above the chosen cutoff is at most . None of these summaries gives conditional correctness for a specific singleton.
Key Parameters
The main parameters that govern the conformal procedure are the miscoverage level , the calibration size , and the choice of nonconformity score.
alpha: The miscoverage level, so is the nominal coverage. It specifies the target marginal miscoverage and determines the interval width through the quantile index. There is no universally correct value; this example uses .n(calibration size): Controls the rank grid and the probability-scale rounding bound. A finite threshold requires ; otherwise the full-space branch is valid. Stability and width also depend on the score distribution.scores(nonconformity score): The main modelling lever. Absolute residuals give constant-width intervals; scaled residuals (dividing by a training-only local noise estimate) give locally adaptive width; signed residuals give one-sided bounds. The same strictly increasing scalar transform applied to every score preserves the sets; an input-dependent normalization generally changes them.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!