can this method work with n_samples < n_features?
how does it compare to ransac which we already have?

thanks

For the n_samples < n_features case it did some modifications in commit 61a5195 in order to fall back to least squares in this case as this Theil-Sen implementation can be seen as a generalization of the ordinary least squares method anyways.
So far I have not found a direct comparison between RANSAC and Theil-Sen but I'm not that familiar with RANSAC that's why I should have deeper look into this subject. So far, my impression is that RANSAC is a more heuristic approach coming from computer science while Theil-Sen comes from robust statistics.

Not to derail the conversation, but is this kind of estimator more appropriate for statsmodels? (My thinking is scikit-learn is focused on predicting outputs while statsmodels is focused on inferences about parameters.)

@chrisjordansquire In order to predict anything you always have to train the parameters of your model and quite often your training samples are less than perfect (i.e. having outliers or if the error is not normally distributed). This is where non-parametric methods like Theil-Sen shine since they do not assume the errors to be normally distributed.
The new RANSAC estimator also seems to fit into the class of robust estimators. From my understanding it takes an estimator like LinearRegression and presents a proper subset to it in order to remove the negative effects of outliers. In this way you could also say that all RANSAC does is to help LinearRegression to infer the right parameters.

Coverage remained the same when pulling b374e7e on FlorianWilhelm:theilsen into 64fd085 on scikit-learn:master.

Not to derail the conversation, but is this kind of estimator more appropriate
for statsmodels?

I am also hesitating on whether this estimator should go in scikit-learn,
or whether we should push users to statsmodels. The reason that we have a
RANSAC is that it is widely used in computer vision, where people really
only about it's ability to fit and predict, and not corresponding
p-values.

I am +0 on including it in scikit-learn: robust models are sometimes
really useful for prediction. However, I would indeed like to hear a
little discussion about the usecases in mind.

@GaelVaroquaux: I understand the hesitation to include yet another estimator especially since RANSAC and Theil-Sen seem to open the gates for a new class of robust estimators and we want scikit-learn to be focused on its goals: An easy-to-use machine learning library with a consistent interface.
But for exactly that reason we are heavily using scikit-learn at Blue Yonder (http://www.blue-yonder.com/en/) to provide predictive analysis products to our customers. At its core everything we do is about making good predictions with the help of machine learning. For our own algorithms we have adapted the scikit-learn interface in order to extend it and that works really well for us.
In some projects we had to deal with suboptimal data (hard to determine outliers) from customers where we applied Theil-Sen successfully. Since Theil-Sen is quite a popular and well-known algorithm my idea was to include it directly into scikit-learn where others can use it too and I got the okay of the management to do so. I am convinced that Theil-Sen (although coming from statistics) is a robust machine learning algorithm that nicely complements Scikit-Learn whenever the data you get is suboptimal in various kinds.

I'm a bit surprised because I never noticed that there is a clean separation between statsmodel and scikit-learn. Clearly this issue here is the wrong place, but I would pretty much appreciate a discussion about this. At least for me, I wouldn't be surprised at all when I would find a Theil-sen implementation in scikit-learn (and another one in scipy and one in statsmodel ;-) )
On the contrary, I heavily use scikit-learn and would be happy to be able to try a Theil-sen (e.g for preprocessing/regularisation) without introducing another dependency (I pretty much never need to introduce statsmodel as a dependency) and with the same interface as the other algorithms.
Last point from me, if someone spends time and wants to contribute to scikit learn, shouldn't the prior be +1 as long as there are no strong arguments against the contribution?
Therefore +1 from me as long as there is no official "that's not scikit-learn" guideline/rule (and that would be great!)

to make your point I would compare performance and speed with the RANSAC we already have. Often computer scientists hacks make things pretty efficient in practice... statisticians care much less about performance ... in general ... :)

Coverage increased (+0.06%) when pulling f92074a on FlorianWilhelm:theilsen into 031a3fc on scikit-learn:master.

@arjoly Regarding the new examples, I used the Unix time command and looked at the user time since the script finishes only when I close the plotting windows.

• plot_theilsen: 1.865s
• plot_robust_fit: 6.152s

I won't be able to find much time for this before next week. @ogrisel and @GaelVaroquaux feel free to finish the review and merged before if you want.

@arjoly @GaelVaroquaux What are the next steps? Do I need to ping someone like @larsmans in order to get this merged? The related PR #3764 has also [MRG+1].

I know nothing about this estimator, so no.

Is this ready to go? @ogrisel and @GaelVaroquaux I see +1's from you, and another implicit +1 by @arjoly !

Both I and @GaelVaroquaux already gave +1 and the last batch of comments by @arjoly have been seemingly addressed, let's merge.

Thanks again @FlorianWilhelm for the contribution (especially the effort in the doc and examples)!

Thanks heaps. This is a quality contribution!

Congratulation @FlorianWilhelm 👍

Thank you all for helping me making my first Scikit-Learn contribution! It was a lot of work but I had tons of fun :-)

@FlorianWilhelm I forgot: can you please open a new PR to add an entry to the doc/whats_new.rst file? Please link your name either to your github account or a personal webpage of your choice at the end of the file.

@ogrisel I updated whats_new.rst in PR #3870.

Stupid question: is there a simple way to make this run fast? It is quite slow in the common tests, even on trivial datasets. I guess setting n_subsamples is the trick, but is only possible if you know the number of samples, right?

@amueller You could also use max_subpopulation to reduce the maximum number of samples that are considered. If you want to reduce the runtime with n_subsamples you need to know the number of samples.

Would it be appropriate to allow n_subsamples as a float to be a proportion of the number of samples?

No, the complexity is $\binom{n_{samples}}{n_{subsamples}}$. If you consider Pascal's triangle, the closer you get to the center of a row the larger the number is. So consider you have n_features and you fit with intercept a problem with n_samples, the efficiency starts to become better only when n_subsamples is larger or equal than n_samples - n_features.
What we could do is to specify it as a ratio of the number of features. If you choose ratio=1.0 (the default) you have the complexity $\binom{n_{samples}}{n_{subsamples}}$. Otherwise we change n_subsample in order to get the complexity $\binom{n_{samples}}{ (n_{samples} - ratio * n_{features}) }$.