[Models] ConformalSeasonalPool: documented n_samples sufficiency rule covers only the lower bound and its example is off by one; no validation or warning when the pool is too small

[shivamlalakiya]

What happened + What you expected to happen

Up front, so this is not mistaken for something bigger than it is: I am not reporting a
coverage-guarantee violation, and I am not reporting the orientation correction itself as
wrong. _oriented_index (models.py:4296-4302) is a genuine finite-sample orientation
correction, its saturation at 0.0/1.0 is a necessary domain constraint (np.quantile rejects
any level outside [0, 1]), the saturated result is still wider than a plain uncorrected
quantile, and the behaviour is both documented and unit-tested. This report is narrower: the
documented sufficiency rule describes only one of the two bounds, its worked example is off by
one, and nothing enforces or warns about the precondition at runtime.

1. The documented rule covers the lower bound only. ConformalSeasonalPool.__init__
docstring, models.py:4268-4270:

n_samples (int, default=100): Number of mixture samples used to estimate prediction intervals.
For a level-L interval, at least ceil(2/(1 - L/100)) - 1 samples are needed before
the orientation-corrected lower bound is non-degenerate (e.g., ≥40 for a 95% interval).

That formula is exactly right — for the lower bound. The upper bound stays pinned to
max(sample) until roughly twice that many samples, and the docstring says nothing about it.
Measured against the installed 2.1.1 (a = 1 - L/100):

level	documented rule ceil(2/a)−1	lower bound non-degenerate at	upper bound non-degenerate at
80	10	10	19
90	20	20	39
95	39	39	79
99	199	199	399

The upper threshold is ceil(4/a) − 1 throughout. So a user who reads the docstring and sets
n_samples=40 for a 95% interval gets a correct lower bound and an upper bound that is still
exactly max(sample).

2. The worked example is off by one. The formula gives ceil(2/0.05) - 1 = 39, and n=39 is
where the lower bound becomes non-degenerate (verified in the table above). The parenthetical
says ≥40.

3. At stock defaults, level=[99] is degenerate on both rails, silently. n_samples
defaults to 100, but level 99 needs 199 for the lower bound and 399 for the upper. So
ConformalSeasonalPool(season_length=12).forecast(y, h, level=[99]) returns literally the min and
max of the 100 mixture draws, and no warning is emitted. The docstring's rule, even if the user
finds it, only tells them they are under 199 — not that the upper bound would need 399.

4. The documented remedy does not apply to predict_in_sample. That path
(models.py:4405-4414) draws from R = model_["calib_residuals"], whose size is fixed by
calib_frac × history — n_samples has no effect on it, so the parameter the docstring points at
cannot fix this path. With 60 observations, season_length=12 and all defaults, R.size is 30,
which is below the 39 needed for level 95, and both offsets saturate to R.min()/R.max().
predict_in_sample's own docstring carries no note about this, and no warning is emitted.

5. n_samples is never validated. __init__ (models.py:4285-4294) validates only
variant. n_samples=1 is reachable from the public API and produces a zero-width interval
(lo == hi exactly).

Expected behaviour

1. The docstring should give both thresholds — ceil(2/a) − 1 for the lower bound and
   ceil(4/a) − 1 for the upper — and the 95% example should read ≥39, matching its own
   formula.
2. The same note should appear on predict_in_sample, phrased against R.size /calib_frac
   rather than n_samples, since n_samples cannot influence that path.
3. __init__ should reject n_samples < 1.
4. Optionally, and most useful to a user who never reads the docstring: warnings.warn when the
   pool size falls below the threshold for a requested level — the formula is already written
   down, so this is a one-line check against a number the code already knows.

I am happy to open the PR for any subset of these. Items 1-3 are small and mechanical; item 4 is
a behaviour addition and therefore your call, which is why this is an issue rather than a PR.

Versions / Dependencies

Click to expand

statsforecast: 2.1.1   (PyPI wheel)
python:        3.13.11
numpy:         2.4.2
pandas:        2.3.3
platform:      macOS-15.7-arm64-arm-64bit-Mach-O

Reproducible example

import math
import warnings

import numpy as np

from statsforecast.models import ConformalSeasonalPool as C

rng = np.random.default_rng(0)
y = (rng.standard_normal(120) + np.tile(np.arange(12), 10)).astype(np.float64)

# (1)+(2) documented rule vs measured thresholds for both bounds
print(f"{'level':>6} {'doc rule':>9} {'lower ok at':>12} {'upper ok at':>12}")
for lv in (80, 90, 95, 99):
    a = 1 - lv / 100
    doc = math.ceil(2 / a) - 1
    lo = max(n for n in range(1, 900) if C._oriented_index(a / 2, n) == 0.0) + 1
    hi = max(n for n in range(1, 900) if C._oriented_index(1 - a / 2, n) == 1.0) + 1
    print(f"{lv:>6} {doc:>9} {lo:>12} {hi:>12}")

# (3) stock defaults, level=99 -> both rails saturated, no warning
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    out = C(season_length=12).forecast(y=y, h=3, level=[99])
    print("\nlevel=99 at defaults, warnings emitted:", [str(x.message) for x in w])
print("oriented lo, hi at n_samples=100, level=99:",
      C._oriented_index(0.005, 100), C._oriented_index(0.995, 100))

# (4) predict_in_sample: pool is R.size, not n_samples
y60 = (rng.standard_normal(60) + np.tile(np.arange(12), 5)).astype(np.float64)
m = C(season_length=12).fit(y60)
R, fitted = m.model_["calib_residuals"], m.model_["fitted"]
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    ins = m.predict_in_sample(level=[95])
    print("\npredict_in_sample warnings emitted:", [str(x.message) for x in w])
lo_off = ins["fitted-lo-95"] - fitted
hi_off = ins["fitted-hi-95"] - fitted
fin = ~np.isnan(lo_off)                      # `fitted` is NaN before calib_start
print(f"R.size={R.size}  lo offset == R.min(): {np.allclose(lo_off[fin], R.min())}"
      f"  hi offset == R.max(): {np.allclose(hi_off[fin], R.max())}")

# (5) n_samples=1 is accepted and gives a zero-width interval
o = C(season_length=12, n_samples=1).forecast(y=y, h=3, level=[95])
print("\nn_samples=1 -> lo == hi:", bool(np.allclose(o["lo-95"], o["hi-95"])))

Output:

 level  doc rule  lower ok at  upper ok at
    80        10           10           19
    90        20           20           39
    95        39           39           79
    99       199          199          399

level=99 at defaults, warnings emitted: []
oriented lo, hi at n_samples=100, level=99: 0.0 1.0

predict_in_sample warnings emitted: []
R.size=30  lo offset == R.min(): True  hi offset == R.max(): True

n_samples=1 -> lo == hi: True

Additional context

I checked for existing reports before filing: _oriented_index and ConformalSeasonalPool
appear only in #1159 (the PR that added them) and #1179 (the CSP tutorial). #1159's own
description documents the orientation correction and the coverage gap it closed, but does not
discuss the saturation thresholds or the missing validation, so this is adjacent prior art rather
than a duplicate. #1119 is a different conformal-interval defect in the ConformalIntervals
wrapper.

---

AI assistance. Very limited, mechanical use of Claude Code to assist with tracing the
quantile and residual call paths in models.py, scaffolding the reproduction script above, and
drafting the boilerplate for this description. The core analysis and logic were human-driven.
Every number quoted is from a run of the snippet as shown. I could not find a stated policy on
AI-assisted contributions in this repository, so I am disclosing by default.

[shivamlalakiya]

One adjacent observation from preparing the reproduction above. Flagging it rather than sitting on it, with the caveat up front that I have not established a direction of harm, so please treat this as a question rather than a defect report.

_oriented_index returns a level, which is then passed to np.quantile. With NumPy's default interpolation (method="linear", Hyndman & Fan type 7) a level q maps to the virtual index q*(n-1), so the level k/n lands at k - k/n rather than on order statistic k:

n	intended k	virtual index hit	order statistic wanted
30	30	29.0000	29
50	49	48.0200	48
100	97	96.0300	96

Recovering order statistic k exactly needs either the level (k-1)/(n-1), or method="higher" on the raw level.

Whether this matters in practice I genuinely do not know. The offset is sub-rank, it shrinks as n grows, and its sign depends on the local spacing of the pooled scores, so it may well be immaterial next to the sampling noise at the pool sizes people actually use. I have not measured the coverage consequence and I am not claiming one.

Raising it only because it is on the same code path as the report above, and because the correction is computed exactly and then delivered through an API whose level-to-index convention is a separate choice. Happy to drop this entirely if it is already understood or considered not worth the churn.

Reproduction: https://doi.org/10.5281/zenodo.21799203 (probes/verify_statsforecast_rebuttal.py, the INTERPOLATION section).

---

AI assistance. Same disclosure as the issue body: limited, mechanical use of Claude Code for tracing call paths and drafting; the analysis is human-driven and every number is from a run of the committed probe.