Conformal prediction intervals for chemists (no PhD required)

Here is a sentence you have seen in papers and slide decks: "the model achieved an R² of 0.91 and a cross-validated RMSE of 0.3 eV." Here is the question nobody asks: what is the error bar on the next prediction?

Not the average error over the last hundred compounds — the error on the compound you are about to screen, the solvent you are about to try, the formulation you are about to mix. An RMSE of 0.3 eV is perfectly consistent with half the predictions being off by 0.1 and the other half by 0.9, and you cannot tell which half the next candidate is in.

Conformal prediction fixes exactly this problem, distribution-free, for any model — random forest, Gaussian process, gradient boosting, even a look-up table. It needs no retraining and about twenty lines of code. This post is the whole method, the coverage check that keeps it honest, and the assumption that chemistry keeps breaking.

The idea in one paragraph

Take a subset of your data the model never trains on. Predict it. Measure how wrong the model was on those points — those residuals are your raw material. Whatever the model was wrong by on the calibration set, it will (with a stated probability) be at most that wrong on the next point. So the interval is: prediction ± the appropriate quantile of the calibration residuals. That is it. No Gaussian assumption, no closed-form derivations, no model-specific machinery.

The formal guarantee: if the calibration points and the new point are exchangeable (same distribution, no ordering), then the interval contains the truth with probability at least 1 − α on average. "On average" is the honest word — some individual intervals will miss, and the long-run miss rate is what is controlled.

The method in twenty lines

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

rng = np.random.default_rng(0)

# Step 0: three-way split. Train / calibrate / (later) test.
X_tr, X_tmp, y_tr, y_tmp = train_test_split(X, y, test_size=0.4, random_state=0)
X_ca, X_te, y_ca, y_te = train_test_split(X_tmp, y_tmp, test_size=0.5, random_state=0)

# Step 1: fit on the training set only.
model = GradientBoostingRegressor(random_state=0).fit(X_tr, y_tr)

# Step 2: nonconformity scores on the calibration set (absolute residuals).
scores = np.abs(y_ca - model.predict(X_ca))

# Step 3: the finite-sample quantile — note the ceil((n+1)(1-alpha)) correction.
alpha = 0.10
n = len(scores)
k = int(np.ceil((n + 1) * (1 - alpha)))
q = np.sort(scores)[min(k, n) - 1]

# Step 4: intervals on new data.
y_hat = model.predict(X_te)
lower, upper = y_hat - q, y_hat + q

Two details that are not optional:

• The `(n+1)` correction. Using the plain 1-alpha quantile of the residuals under-covers at small n. With the correction, a 90% interval on a 20-point calibration set uses the second-worst residual — wide, honest, and nominally valid.
• Never calibrate on training residuals. Training residues are optimistically small; calibrated intervals would under-cover silently. The calibration set must be held out from fitting.

The check that keeps it honest

Run the interval on a test set and count the coverage. Do this every time you retune anything:

covered = ((y_te >= lower) & (y_te <= upper)).mean()
width = (upper - lower).mean()
print(f"target coverage {1 - alpha:.2f} | observed {covered:.2f} | mean width {width:.3f}")

If observed coverage is well below nominal on data drawn from the same distribution, something is wrong in the split or the quantile. If it is much higher than nominal (say 0.99 for a 0.90 target), your intervals are wider than they need to be — not a correctness failure, but a loss of resolution. Coverage and width are the two numbers that matter together; a method that games either one is a method that lies.

Where chemistry breaks the guarantee

Exchangeability is the assumption, and experimental science is where it dies:

1. New chemistry is not exchangeable with old chemistry. A calibrationset of small molecules does not cover a novel scaffold. Coverage degrades exactly where the prediction is most valuable. This is not fixable with math — it is a data problem. The honest response is domain checks and flags (see below), and treating intervals as lower bounds on uncertainty outside the domain.

1. Batch effects and drift. Solvent lots, instrument recalibration, seasonal humidity. If the calibration set is old, the residuals are not representative of today's measurement process.

1. Heteroscedasticity. A single ± q gives every prediction the same width. If your model is much more uncertain on some compounds, use the normalized score: divide each calibration residual by a model-provided uncertainty (a Gaussian process's predictive std, or a k-NN distance to the training set), take the quantile of those ratios, and scale each interval by its own uncertainty estimate. Same code, two extra lines, and intervals that breathe.

1. Tiny calibration sets. With n = 10, the 90% interval is set by the worst residual — a single outlier widens everything. That is the correct behaviour, and it is also a message: go measure more calibration data before making decisions. Splitting a 30-row dataset into train/calibrate/test leaves almost nothing; prefer cross-conformal variants (calibrate across CV folds) when data is this thin.

What this looks like in practice

Matflow's prediction workflow uses split-conformal intervals on every prediction, alongside an extrapolation flag that marks outputs outside the training envelope — the engineering response to the exchangeability problem above. Benchmark reports and their reproduction bundles are public at /benchmarks, and the protocol is documented at /benchmarks/methodology. If you are running your own stack, everything in this post is scikit-learn and numpy, and the small-data tutorial covers the evaluation discipline that goes with it.

The reframe to take away: the interval is the product, not the point estimate. A prediction of 3.2 eV means nothing; "3.2 eV, 90% interval [2.9, 3.5], in-domain" is a sentence you can make a decision with. When the data does not support a sharp answer, the interval says so — and a model that admits uncertainty is worth more than one that hides it.

Honest limits

• Conformal intervals are marginal, not conditional: they are right on average, but a specific subgroup may be systematically under-covered. Check coverage per material family or process window when the stakes justify it.
• The guarantee is only as good as the exchangeability assumption; domain shift silently erodes it. That is why the flag exists and why "in-domain" belongs next to every interval.
• Conformal prediction does not improve accuracy. It quantifies the accuracy you have. Pair it with honest validation, not instead of it.
• Adaptive intervals inherit the quality of the model's own uncertainty estimate — a badly calibrated σ makes badly scaled intervals even if marginal coverage holds.

If you want to go deeper on the evaluation side — what to measure, and what the numbers do and do not support — start with /benchmarks/methodology and the conformal method note inside the prediction docs.