fix: use the finite-sample conformal quantile for conformal intervals

[Nicoló Angileri (nicoloangileri)]

Summary

ConformalIntervals takes plain empirical quantiles of the n = n_windows conformity scores, so its intervals under-cover for any number of windows:

• conformal_error uses np.quantile(cs, lv / 100), a linear interpolation between scores.
• conformal_distribution uses the alpha / 2 and 1 - alpha / 2 quantiles of the 2n pooled values mean -/+ cs.

Split conformal prediction needs the finite-sample correction: the ceil((n + 1) * q)-th smallest score for an upper quantile and the floor((n + 1) * q)-th for a lower one. ConformalSeasonalPool already applies it (_oriented_index), the shared ConformalIntervals path does not. With 10 windows, a 90% interval covered 82% before this change and covers 91% after it.

Change

• New _conformal_quantile in utils.py, next to ConformalIntervals. It returns the exact order statistic of _conformal_rank(n, q, upper): ceil((n + 1) q) for an upper bound and floor((n + 1) q) for a lower one, clipped to [1, n]. The tail is passed explicitly rather than inferred from q, because an absolute-score margin is an upper bound even at levels below 50. There is no interpolation, because interpolating at floor((n + 1) q) / n lands slightly above the target score and under-covers the lower tail.
• conformal_error and conformal_distribution call it instead of np.quantile. conformal_error always passes upper=True. conformal_distribution keeps its original alpha / 2 and 1 - alpha / 2 cuts on the pooled scores, with upper=False for the lower cuts and upper=True for the upper ones.
• ConformalSeasonalPool calls the same helper in _intervals_from_samples and predict_in_sample, and _oriented_index is removed. Its lower bounds now use the exact order statistic, and its upper bounds use the exact ceil((n + 1) q)-th sample instead of a point slightly above it.

With absolute conformity scores, both ConformalIntervals methods now return mean -/+ the same order statistic; for conformal_distribution it is occasionally one rank higher, since the pooled 2n values use a slightly more conservative rank.

Known limitation, not addressed here

The default n_windows=2 cannot reach typical levels. With two windows no finite interval covers more than 2/3, so the 80/90/95 bounds all use the largest score. That is the correct behavior for two scores, but the default itself stays too small. It would be better handled separately, for example with a warning or a larger default.

Coverage

Setup: 4000 i.i.d. N(0, 1) series of length 100, HistoricAverage with ConformalIntervals(h=1, n_windows=k). Coverage is measured on the next observation, and the Monte Carlo standard error is about 0.005. Both methods give the same numbers after the change.

Two-sided coverage, conformal_distribution:

n_windows	80% before	80% after	90% before	90% after	95% before	95% after
2	0.598	0.672	0.638	0.672	0.655	0.672
5	0.689	0.830	0.774	0.830	0.810	0.830
10	0.736	0.814	0.823	0.906	0.875	0.906
20	0.765	0.804	0.855	0.904	0.905	0.950
50	0.781	0.801	0.879	0.897	0.928	0.961

Each tail after the change (share of observations below lo / above hi, target (1 - level) / 2):

n_windows	80%: lo / hi	90%: lo / hi	95%: lo / hi
10	0.094 / 0.092	0.044 / 0.050	0.044 / 0.050
20	0.100 / 0.096	0.046 / 0.051	0.025 / 0.025
50	0.097 / 0.102	0.051 / 0.052	0.017 / 0.022

Levels below 50 with 10 windows (both methods, after the change; the target is the ceil rank ceil(11 * lv / 100) / 11):

level	main (conformal_error)	after	ceil-rank bound
30	0.336	0.364	0.364
40	0.419	0.458	0.455
49	0.493	0.546	0.545

conformal_error before the change, at 90%: 0.649 / 0.782 / 0.831 / 0.859 / 0.880 for 2 / 5 / 10 / 20 / 50 windows.

Tests

• tests/test_utils.py: _conformal_rank for both tails and the clipping; _conformal_quantile returns exact order statistics, handles sequences of levels and other axes, and rejects empty scores; Monte Carlo checks that each tail misses at most (1 - level) / 2 for 10, 20 and 50 scores, and that an upper margin at 30, 40 and 49 covers at least its level (within 3 standard errors).
• tests/test_models.py: with 10 windows, both methods put the bounds exactly at -/+ the 8th and 9th smallest scores for 70% and 80%, and at the 4th, 5th and 6th for 30%, 40% and 49%; the 2-window exact-value check is unchanged.
• tests/test_csp.py: the orientation tests now check the exact order statistics through the shared helper.
• Results: tests/test_utils.py, tests/test_csp.py and tests/test_models.py 261 passed, tests/test_core.py 60 passed. ruff check 0.2.1 passes, and mypy reports no new errors.

Merged main into the branch.

[CLAassistant]

All committers have signed the CLA.

[Nicoló Angileri (nicoloangileri)]

Hi Saul (@nasaul), a gentle ping on this PR when you get a chance. The CLA is signed, and I'm happy to change the approach if you'd prefer a different one. Thanks!

[Nicoló Angileri (nicoloangileri)]

Hi Saul (@nasaul), a gentle follow-up on this one (and its twin in neuralforecast, Nixtla/neuralforecast#1628). Happy to adjust the approach if you'd prefer something different. Thanks!

[Saul (nasaul)]

Thanks for this, Nicoló Angileri (@nicoloangileri), and sorry for the late reply. The diagnosis is right and the fix is the right direction.

One request on the motivation: the description mixes two separate problems, and this PR solves only one of them.

1. The quantile isn't finite-sample corrected. The plain lv / 100 quantile under-covers for any n_windows. This PR fixes that.
2. The default n_windows=2 is too small to reach typical levels. With two windows, no finite interval covers more than 2/3, so the 80/90/95 bounds all collapse to the largest score. That's correct behavior, but the problem itself remains.

Leading with "a 90% interval covers about 64% with the default n_windows=2" suggests this PR fixes the default case, and it can't. Could you reframe the motivation around (1), and mention (2) as a known limitation to handle separately?

On the implementation: the statistical correction should be the one ConformalSeasonalPool already applies (_oriented_index), not a second variant. Could you move that correction into a shared helper in utils.py (next to ConformalIntervals) and have both ConformalIntervals methods and CSP call it? That way there's one implementation of the finite-sample quantile, one set of tests, and the same guarantee everywhere.

While extracting the helper, please check the lower tail. With linear interpolation, floor((n + 1) q) / n lands slightly above the j-th residual, and in simulation that tail under-covers (about 17% misses against a 10% target at n=10 / 80%). Taking the exact order statistic would fix it, using either an index into the sorted scores or method="inverted_cdf".

A rebase on main will also be needed. Thanks again for the thorough write-up and coverage tables!

[Nicoló Angileri (nicoloangileri)]

Thanks a lot for the careful review, Saul (@nasaul). All points are addressed in a18c71d:

1. Motivation: the description now focuses on the finite-sample correction, and the n_windows=2 default is listed as a known limitation to handle separately.
2. Single implementation: the correction lives in _conformal_quantile in utils.py, next to ConformalIntervals. Both ConformalIntervals methods and ConformalSeasonalPool (_intervals_from_samples and predict_in_sample) call it, and _oriented_index is gone. conformal_distribution keeps its original pooled cuts and only swaps np.quantile for the helper.
3. Lower tail: the helper takes the exact order statistic, with no interpolation. In the simulation, each tail now misses close to its target: at 80% with 10 / 20 / 50 windows the lower tail misses 0.094 / 0.100 / 0.097. There is also a Monte Carlo test for both tails in tests/test_utils.py.
4. Main: I merged main into the branch instead of rebasing, to avoid a force-push. Happy to rebase if you prefer a linear history.

One side effect worth flagging: with absolute conformity scores, both ConformalIntervals methods now return the same bounds (the pooled one is occasionally one rank higher), and CSP's upper bounds move to the exact order statistic, slightly below the interpolated point they used before. The coverage tables in the description are updated.

[Saul (nasaul)]

Great job on the update, Nicoló Angileri (@nicoloangileri)! The reframed motivation reads well, the shared helper is exactly what we were after, and the exact order statistic fixes the CSP lower tail. Coverage now sits on the min(k, n)/(n + 1) bound in our simulations.

One last thing before merging: conformal_error now under-covers for levels below 50.

_conformal_rank infers the tail from q: floor when q < 0.5, ceil otherwise. That's right for CSP, whose lower and upper quantiles are always on opposite sides of 0.5. But conformal_error passes lv / 100 for an absolute-score margin, which is always an upper quantile. For lv < 50 it takes the floor path and lands one rank too low. With 10 windows:

level	main	this PR	ceil rank
30	0.336	0.272	0.362
40	0.419	0.362	0.454
49	0.493	0.454	0.544

At level 40, floor(11 * 0.4) = 4 picks the 4th smallest score, which covers ~4/11; the ceil rank of 5 is needed. This is a regression vs main, which met the target at those levels.

Could you make the tail explicit instead of inferring it from q? For example _conformal_rank(n, q, upper: bool), with conformal_error always passing upper=True and CSP passing upper=False / True for its lower and upper bounds. Please also check that conformal_distribution gets the upper-tail rank for low levels, and add a test with a level below 50.

Thanks again, this is a really nice piece of work.

[Nicoló Angileri (nicoloangileri)]

Good catch, thank you! Inferring the tail from q was wrong for absolute-score margins. Fixed in be69272:

• _conformal_rank(n, q, upper) and _conformal_quantile(scores, q, upper, axis) now take the tail explicitly (a bool, or one bool per level).
• conformal_error always passes upper=True. conformal_distribution passes upper=False for its alpha / 2 cuts and upper=True for its 1 - alpha / 2 cuts; those cuts stay on their own side of 0.5 at every level, so low levels already got the upper rank there, and now it is explicit. CSP passes False / True for its lower and upper bounds.
• New tests for levels 30, 40 and 49 with 10 windows: both methods return exactly -/+ the 4th, 5th and 6th smallest scores, plus a Monte Carlo coverage check for an upper margin below 50.

Same setup as your table (10 windows), coverage after the fix for both methods: 0.364 / 0.458 / 0.546 at 30 / 40 / 49, on the ceil-rank bound. The description is updated with these numbers.

[Saul (nasaul)]

Looks good, just waiting on the CI to run

[Saul (nasaul)]

LGTM, thanks for the contribution, happy forecasting