DEV Community

Devanshu Biswas
Devanshu Biswas

Posted on

The Theil-Sen estimator: fit a robust regression line with nothing but a median

Ordinary least squares picks the slope that minimises the sum of squared residuals. That squaring is its fatal flaw: a single far-off point has an enormous squared error, so it pivots the whole line toward itself. OLS's breakdown point is zero — one bad point is enough. The Theil-Sen estimator answers with the humblest robust idea in statistics, and it fits a line using nothing more exotic than the median.

The whole estimator, in two medians

Two points define a line, and its slope is rise over run. So each pair of your data points is a tiny, independent slope estimate. Theil-Sen computes the slope of every pair and takes the median. The intercept is then just median(yᵢ − m·xᵢ).

function pairSlope(p, q) {
  const dx = q.x - p.x;
  if (Math.abs(dx) < 1e-9) return null;   // same x -> undefined slope, skip
  return (q.y - p.y) / dx;                // rise over run
}

function allPairSlopes(pts) {
  const slopes = [];
  for (let i = 0; i < pts.length; i++)
    for (let j = i + 1; j < pts.length; j++) {   // each unordered pair once
      const s = pairSlope(pts[i], pts[j]);
      if (s !== null) slopes.push(s);
    }
  return slopes;                                  // n(n-1)/2 of them
}

function theilSen(pts) {
  const m = median(allPairSlopes(pts));            // slope = median of slopes
  const intercepts = pts.map(p => p.y - m * p.x);  // b that centres each point
  const b = median(intercepts);                    // intercept = their median
  return { m, b };                                 // y = m*x + b, done
}
Enter fullscreen mode Exit fullscreen mode

That's it. No loss function, no iteration, no gradient, no learning rate, no threshold to tune, and no random sampling. The estimate is exact and deterministic — run it twice on the same data and get the same line to the last bit.

Why the median beats OLS on outliers

The median depends only on the rank of a value, not its magnitude. You can send half-minus-one of the numbers to infinity and it doesn't move.

Now think about what an outlier does to the pairwise slopes. It can only corrupt the pairs it belongs to — a minority of all pairs — and those corrupted pairs produce extreme slopes that land in the tails of the distribution. The mean (which OLS-style averaging chases) gets dragged toward the tail. The median only cares which side of centre a value sits on, so it steps right over the tail slopes.

The ~29% breakdown point

How much contamination can Theil-Sen take? The median of a list breaks at 50%. But Theil-Sen's slope is a median of pairs, and a pair is clean only if both its points are inliers. If a fraction ε of points are bad:

P(clean pair) = (1 - ε)²
median tips once bad pairs reach half:
  (1 - ε)² = 1/2   ->   ε = 1 - 1/sqrt(2) = 0.2929...
Enter fullscreen mode Exit fullscreen mode

So Theil-Sen tolerates just under 29.3% arbitrary contamination — a colossal jump from OLS's 0% — while staying a simple closed formula. (RANSAC and least-median-of-squares can push past 50%, but they need a threshold and random sampling; Theil-Sen needs neither.)

Where it sits, and where it still bends

Its one cost is the O(n²) sweep over pairs, trivial for the small, low-dimensional trend problems where it shines (it's a favourite for climate and econometric trend estimation). In practice you reach for the library versions:

# SciPy — slope, intercept, and a slope confidence interval
from scipy.stats import theilslopes
slope, intercept, lo, hi = theilslopes(y, x, 0.95)

# scikit-learn — multivariate Theil-Sen (subsamples in high-D)
from sklearn.linear_model import TheilSenRegressor
ts = TheilSenRegressor(random_state=0).fit(X, y)
Enter fullscreen mode Exit fullscreen mode

Theil-Sen is the elegant middle ground between the smooth all-points estimators (OLS, Huber) and the sample-and-vote world of RANSAC: more robust and simpler to use than Huber, more deterministic and tuning-free than RANSAC. Its one blind spot — shared with almost every method here — is a high-leverage point far out in x, which can keep a small residual and still tilt the pairwise slopes.

Reach for it when your data is one- or low-dimensional, you expect up to a quarter of it to be junk, and you want a robust slope with literally nothing to tune.

Drag points, add outliers, and watch the median-of-slopes line hold while least-squares gets dragged away: https://dev48v.infy.uk/ml/day53-theil-sen.html

Top comments (0)