BUG: CovDetMCD/CovDetS/CovDetMM raise ValueError when nobs < 2*k_vars + 4
Describe the bug
CovDetMCD, CovDetS and CovDetMM raise ValueError: Percentiles must be in the range [0, 100] whenever the data satisfy nobs < 2 * k_vars + 4, regardless of the arguments passed.
The failure comes from the deterministic starting-set construction, not from the estimator itself, so it happens before any user parameter is consulted. In statsmodels/robust/covariance.py, _cov_starting (line 1398 on main):
percentiles = [(k_vars + 2) / nobs * 100 * 2, 25, 50, 85]
cutoffs = np.percentile(d, percentiles)The first element is 200 * (k_vars + 2) / nobs, which exceeds 100 — and so trips numpy's range check — exactly when nobs < 2 * k_vars + 4.
The affected region is inside the documented usable range of the estimators. CovDetMCD requires only h > k_vars for the subset covariance to be non-singular, and h defaults to well above k_vars here; e.g. at nobs=60, k_vars=30 the max-breakdown h = (nobs + k_vars + 1) // 2 = 45 > 30, so the estimate is well defined. R's robustbase::covMcd(x, nsamp = "deterministic") computes the same case successfully, emitting only a n < 2 * p warning.
Code Sample
import numpy as np
from statsmodels.robust.covariance import CovDetMCD, CovDetS, CovDetMM
x = np.random.default_rng(0).standard_normal((60, 30)) # nobs < 2*k_vars + 4 = 64
CovDetMCD(x).fit(45) # ValueError: Percentiles must be in the range [0, 100]
CovDetS(x).fit() # same
CovDetMM(x).fit() # same
x2 = np.random.default_rng(0).standard_normal((60, 28)) # nobs >= 2*k_vars + 4 = 60
CovDetMCD(x2).fit(44) # OK
CovDetS(x2).fit() # OK
CovDetMM(x2).fit() # OKTraceback:
File "statsmodels/robust/covariance.py", line 2189, in fit
starts = _get_detcov_startidx(z, h_start, options_start)
File "statsmodels/robust/covariance.py", line 1618, in _get_detcov_startidx
cov_all = _cov_starting(z, standardize=False, quantile=0.5)
File "statsmodels/robust/covariance.py", line 1399, in _cov_starting
cutoffs = np.percentile(d, percentiles)
ValueError: Percentiles must be in the range [0, 100]The boundary is exactly as the formula predicts. Sweeping k_vars for several nobs and recording the first value that raises:
nobs |
predicted first failing k_vars (> (nobs-4)/2) |
observed |
|---|---|---|
| 41 | 19 | 19 |
| 60 | 29 | 29 |
| 61 | 29 | 29 |
| 100 | 49 | 49 |
| 200 | 99 | 99 |
| 301 | 149 | 149 |
It is not a matter of how fit is called. At nobs=60, k_vars=30, all of h ∈ {31, 35, 40, 45, 46, 50, 55, 59} raise identically, as do reweight=False, h_start ∈ {45, 50, 60} and maxiter_step=1 — consistent with the expression depending only on nobs and k_vars.
Expected Output
An estimate, as in the k_vars=28 control above and as robustbase produces for the same input.
The intent of the first percentile appears to be "take roughly the 2*(k_vars + 2) observations closest to the centre as one candidate starting subset". When that count meets or exceeds nobs, the sensible reading is all observations rather than an error, so clipping would preserve the intent:
percentiles = [min((k_vars + 2) / nobs * 100 * 2, 100), 25, 50, 85]I have not opened a PR since I do not know whether the other three percentiles are meant to be adjusted in that regime as well, or whether the maintainers would prefer an explicit guard with a clearer message. Happy to submit one if the clip is the preferred fix.
Context
Found while using CovDetMCD on narrow column subsets of a wide design (a leverage-detection scheme that fits many nobs x p_star submatrices), where p_star near nobs/2 is a normal operating point. Falling back to robustbase was the workaround.
Python: 3.12.13.final.0
OS: Darwin 25.6.0 arm64
statsmodels: 0.15.0
numpy: 2.5.2
scipy: 1.18.1
pandas: 3.0.5
patsy: 1.0.3Line numbers above are from main as of this writing; covariance.py:1398 still contains the unclipped expression.
Source: statsmodels/statsmodels