We're using the Theil-Sen estimator for fitting a line (orange) through mostly straight data. It works very well compared to the linear trendline (blue) in Excel:


We also need to fit mostly curved data with outliers. Is there something similar to the Theil-Sen estimator we could use? Excel's second degree polynomial trendline (blue) for all points, results in a bad fit due to the outliers on the right. When the data is manually truncated to 0 < x < 1.9, the polynomial trendline (orange) fits a lot better. But the outliers may be present in the middle of the data as well.


A second degree polynomial fits our data good enough, but other curves may be fitted as well. Our goal is to get the y position and slope of the fitted curve at each x.

Excel data available here.

  • $\begingroup$ Why not a linear spline (for the seemingly linear data)? $\endgroup$
    – Firebug
    Commented Dec 8, 2017 at 12:46
  • 5
    $\begingroup$ As I understand it, this estimator fits every combination of two data points to first order polynomial - a straight line - and then uses a median of the coefficient values as the result. If so, then a second order polynomial equivalent would fit every combination of three points to a parabola and then find the median of those coefficient values. I have never tried this, though it seems plausible as a curved analog to the standard first order regression. $\endgroup$ Commented Dec 9, 2017 at 20:33
  • $\begingroup$ That sounds quite promising, I will try that tomorrow. $\endgroup$
    – Anlo
    Commented Dec 10, 2017 at 8:00
  • 1
    $\begingroup$ I have both of Theil's original articles. They are on polynomial regression and not linear regression. They are a pill to read, but there is no restriction on how to do it. $\endgroup$ Commented Dec 18, 2017 at 19:00
  • $\begingroup$ @James I believe methods like this have been applied. As a practical matter, I would consider choosing subsets of three points that are spread across the entire dataset and eliminate those that have short ranges or include pairs of nearby points. This is more of an issue for polynomials than linear functions, because a much larger proportion of all triples will be "local" in this sense, whereas you want a global estimator; and it would make the calculation much more efficient. $\endgroup$
    – whuber
    Commented Mar 25 at 16:39

2 Answers 2


James Phillips' suggestion on how to expand the Theil-Sen algorithm to a second degree polynomial worked surprisingly well. There were 762 (x,y)-points in the dataset I tested. Selecting three different points from the 762 can be made in 73 million ways, so instead I put the points into groups of 11 and calculated the median x och y values of each group. $\binom{762/11}{3} = \binom{69}{3} \approx 52000$ which is a more reasonable number of combinations to use.

For each combination of three points, I calculated the $a$ coefficient for $y = ax^2 + bx + c$ using:

$a = \frac{x_3(y_2 - y_1) + x_2(y_1 - y_3) + x_1(y_3 - y_2)}{(x_1 - x_2)(x_1 - x_3)(x_2 - x_3)}$

Using the median of all $a$'s, for each combination of two points I calculated the slope of the line $y - ax^2 = bx + c$ using:

$b = \frac{(y_2 - ax_2^2) - (y_1 - ax_1^2)}{x_2 - x_1}$

Finally, using the the medians of all $a$'s and $b$'s, for each point I calculated the intercept of $y - ax^2 - bx = c$ using:

$c = y_n - ax_n^2 - bx_n$

I used the Remedian median value approximation and the algorithm implemented in C# took ~15 ms to run for 69 points. If the initial median filter is reduced to 5 points, the execution time increases to ~110 ms which still is ok. Running the algorithm on all 762 points results in an execution time of 13 seconds.

Even with the fast 11 point filter, the results looks very good:

2nd degree Theil-Sen

  • 1
    $\begingroup$ At the risk of using particularly arcane technical jargon, "Totally Awesome - for sure, dude!" - a most excellent outcome. $\endgroup$ Commented Dec 29, 2017 at 17:21

According to Dang et al., you should probably take the spatial median of the coefficients instead of the median along each axis. The two are not exactly equivalent. The former minimizes the sum of $L_2$ distances and latter the sum of $L_1$ distances.

This is also the approach taken by the TheilSenRegressor in scikit-learn (see here).

Xin Dang et al., "Theil-Sen Estimators in a Multiple Linear Regression Model", 2009.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.