Repository: scikit-learn/scikit-learn
Pull request: FIX adjust inliner criterion in RANSACRegressor
#### Reference Issues/PRs
Fixes #19497. 
Fixes #16147.

#### What does this implement/fix? Explain your changes.
This fix implement solution discussed in #16147: 
inlier_mask_subset = residuals_subset <= residual_threshold
that replace the current <

This make the code work also with perfect horizontal lines (see #19497). 

Basically I suggest to change the definition of what is an inlier respectivly to residual_threshold. Exactly on the threshold, I suggest to say it is an inlier and not an outliner.

This make a couple of "normal cases" work: horizontal line #19497 (currently we say all the point on a line are outlier!) and other examples described in #16147

Unfortunately This patch make a test fail, but we can keep it with the following modification. The test construct an example where there is no inlier. It is based on the fact when residual_threshold =0 an exception is always thrown. To do this in a clean way we should construct here a real example where there is no inlier. Such an example is not easy to find. 
I suggest to construct a more simple test as before : set the residual_threshold to an other corner value that make always the comparison fail : a nan

I also implement the test described in #19497 (horizontal lines)
#### Any other comments?


Problem statement:
RANSAC - residual threshold calculation
Within the RANSAC algorithm, a residual threshold is calculated:

`residual_threshold = np.median(np.abs(y - np.median(y))`

If more than half of the values of y are equal to the median of y, this returns a residual threshold of 0.  In that case, the line 

`inlier_mask_subset = residuals_subset < residual_threshold`

always returns zero inliers, causing a value error since `inlier_mask_best` is always None.

RANSAC - perfect horizontal line generate an exception
#### Describe the bug
points sampled from an horizontal line without noise generate an exception 

#### Steps/Code to Reproduce
```python
    import numpy as np
    from sklearn.linear_model import LinearRegression, RANSACRegressor
    X = np.arange(100)[:, None]
    y = np.zeros((100, ))

    base_estimator = LinearRegression()
    ransac_estimator = RANSACRegressor(base_estimator, random_state=0)
    ransac_estimator.fit(X, y)

    assert_allclose(ransac_estimator.estimator_.coef_, 0.)
    assert_allclose(ransac_estimator.estimator_.intercept_, 0.)
```

#### Expected Results
No error is thrown

#### Actual Results

```python-traceback
Traceback (most recent call last):
  File "<ipython-input-2-84e926dec73f>", line 8, in <module>
    ransac_estimator.fit(X, y)
  File "/Users/ogrisel/code/scikit-learn/sklearn/linear_model/_ransac.py", line 439, in fit
    raise ValueError(
ValueError: RANSAC could not find a valid consensus set. All `max_trials` iterations were skipped because each randomly chosen sub-sample failed the passing criteria. See estimator attributes for diagnostics (n_skips*).
```

#### Versions
In [1]: import sklearn; sklearn.show_versions()

System:
    python: 3.8.5 (default, Sep  4 2020, 02:22:02)  [Clang 10.0.0 ]
executable: /usr/local/Caskroom/miniconda/base/envs/sklearndev/bin/python
   machine: macOS-10.14.6-x86_64-i386-64bit

Python dependencies:
          pip: 21.0.1
   setuptools: 52.0.0.post20210125
      sklearn: 1.0.dev0
        numpy: 1.19.2
        scipy: 1.6.0
       Cython: 0.29.21
       pandas: None
   matplotlib: 3.3.4
       joblib: 1.0.1
threadpoolctl: 2.1.0

Built with OpenMP: False

Hints from issue discussion:
```python
np.median(np.abs(y - np.median(y))
```
is the definition of the Median Absolute Deviation. For the case that you provide the MAD is equal to zero and this is correct.

So I think that we could do several things to improve the estimator:

1. raise a proper error message when we compute the MAD and explain that one should use another robust (maybe less robust) variability estimator.
2. allow `residual_threshold` to be a callable instead of expecting a floating number.
It's me again @glemaitre. This issue looks interesting. May I work on it?
Yes go for it ;)

Sent from my phone - sorry to be brief and potential misspell.



Thanks. PR soon
The error happens when the RANSACRegressor is created without passing a residual_threshold. In this case RANSAC needs to calculate the residual_threshold from the y vector.
The error thrown during fit is:

```
File "/Users/ashacked/dev/python/scikit-learn/sklearn/linear_model/_ransac.py", line 431, in fit

"RANSAC could not find a valid consensus set. All"

ValueError: RANSAC could not find a valid consensus set. All `max_trials` iterations were skipped because each randomly chosen sub-sample failed the passing criteria. See estimator attributes for diagnostics (n_skips*).
```

The same error will be thrown if we create the RANSACRegressor with  residual_threshold=0.0, and this behaviour is tested in **test_ransac.test_ransac_resid_thresh_no_inliers()**.

And actually passing a Y vector where most of the values are the same (and not passing residual_threshold=0.0) has the same signifficance because then the residual_threshold will be computed from the Y vector and will be again zero. 
So in short the error generated by this scenario is thrown because the input creates a residual_threshold equal to zero and the current implementation does not tolerate that, as the regression test **test_ransac_resid_thresh_no_inliers()** clearly demonstrates.

So one possible resolution could be to accept the current behaviour which seems to be "as designed"

On the other hand, we could ask this question. Should we refuse to handle the case where residual_threshold=0.0. This is actually the case when the input to the Regressor is an horizontal line. gradient is zero.

Having more than 50 percent of target values equal to the median, is actualy a reasonable (although trivial) case. Trying to find a regression model for an horizontal line, will pass a Y vector where all  Y values are the same. There are no outliers here, the gradient is 0, but maybe the RANSAC model should be able to cope with the trivial border case.

The line in [RANSACRegressor.fit](http://ransacregressor.fit)() that fails the case where residual_threshold is 0, is:
`inlier_mask_subset = residuals_subset < residual_threshold`
if we had a < = instead of =, fit() would have worked. Meaning accepting a sample where the predict target value is equal to the sample target value. Which does not seem far fetched to me.
Actually in the Wikipedia page for RANSAC

[https://en.wikipedia.org/wiki/Random_sample_consensus](https://en.wikipedia.org/wiki/Random_sample_consensus), the MATLAB implementation has
 a " < = " and not a " < "

`inlierIdx = find(abs(distance)<=threshDist) 
`

I made this change and worked with 
`inlier_mask_subset = residuals_subset < = residual_threshold`
Ran test_ransac.py and saw that all tests passed except test_ransac_resid_thresh_no_inliers()  which failed as expected because it specifically verifies that we dont accept residual_threshold=0.0 and because of my change we did.
@glemaitre , @jnothman  what do you think?
@alexshacked I made the same change (from "<" to "<=" in my local copy and have found it to work well in practice.  

I would argue that this should be the default behavior.  While trivial, it is a case that can be observed in practice.  I encountered this while looping through client accounts in which the target quantity is constant (or near constant) for a small fraction of accounts.  When the target quantity is constant, it's easy to throw out the case, but it may still be of interest to run the regression when 51% of the observations have the same target value.

Commit message: Merge branch 'main' of https://github.com/scikit-learn/scikit-learn into bugfix_ransac_residual_threshold

Patch to review:
diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst
index f1f376dc641c9..7ed8cb32d8973 100644
--- a/doc/modules/linear_model.rst
+++ b/doc/modules/linear_model.rst
@@ -1265,8 +1265,8 @@ Each iteration performs the following steps:
    whether the estimated model is valid (see ``is_model_valid``).
 3. Classify all data as inliers or outliers by calculating the residuals
    to the estimated model (``base_estimator.predict(X) - y``) - all data
-   samples with absolute residuals smaller than the ``residual_threshold``
-   are considered as inliers.
+   samples with absolute residuals smaller or equal than the 
+   ``residual_threshold`` are considered as inliers.
 4. Save fitted model as best model if number of inlier samples is
    maximal. In case the current estimated model has the same number of
    inliers, it is only considered as the best model if it has better score.
diff --git a/doc/whats_new/v1.0.rst b/doc/whats_new/v1.0.rst
index 602d4b1246878..d4f07988f1b34 100644
--- a/doc/whats_new/v1.0.rst
+++ b/doc/whats_new/v1.0.rst
@@ -217,6 +217,19 @@ Changelog
   :pr:`19426` by :user:`Alexandre Gramfort <agramfort>` and
   :user:`Maria Telenczuk <maikia>`.
 
+- |Fix| :class:`RANSACRegressor` Points with a residual equal to 
+  ``residual_threshold`` member are considered as an inlier.
+  Motivation for this change: We want a perfect dataset to be considered as 
+  inliers, specifically with default parameters or ``residual_threshold`` 
+  equal to zero.
+  Previously the code launched an exceptions with perfect model.
+  Such a comportement also happen with MAD and half of the 
+  points perfect.
+  With this definition, ``residual_threshold`` set to zero no longer set all 
+  points to outliers. Now you can do that by setting ``residual_threshold`` to 
+  nan.
+  :pr:`19499` by :user:`Gregory Strubel <gregorystrubel>`.
+
 - |Efficiency| The implementation of `fit` for `PolynomialFeatures` transformer
   is now faster. This is especially noticeable on large sparse input.
   :pr:`19734` by :user:`Fred Robinson <frrad>`.
diff --git a/sklearn/linear_model/_ransac.py b/sklearn/linear_model/_ransac.py
index 2fc8143f432c8..409edb3fa0b1f 100644
--- a/sklearn/linear_model/_ransac.py
+++ b/sklearn/linear_model/_ransac.py
@@ -95,7 +95,8 @@ class RANSACRegressor(MetaEstimatorMixin, RegressorMixin,
     residual_threshold : float, default=None
         Maximum residual for a data sample to be classified as an inlier.
         By default the threshold is chosen as the MAD (median absolute
-        deviation) of the target values `y`.
+        deviation) of the target values `y`. Residual strictly equal to the
+        threshold is considered valid.
 
     is_data_valid : callable, default=None
         This function is called with the randomly selected data before the
@@ -393,7 +394,7 @@ def fit(self, X, y, sample_weight=None):
             residuals_subset = loss_function(y, y_pred)
 
             # classify data into inliers and outliers
-            inlier_mask_subset = residuals_subset < residual_threshold
+            inlier_mask_subset = residuals_subset <= residual_threshold
             n_inliers_subset = np.sum(inlier_mask_subset)
 
             # less inliers -> skip current random sample
diff --git a/sklearn/linear_model/tests/test_ransac.py b/sklearn/linear_model/tests/test_ransac.py
index 857696bf387d5..6a58e7af92626 100644
--- a/sklearn/linear_model/tests/test_ransac.py
+++ b/sklearn/linear_model/tests/test_ransac.py
@@ -151,12 +151,12 @@ def test_ransac_predict():
 
 
 def test_ransac_resid_thresh_no_inliers():
-    # When residual_threshold=0.0 there are no inliers and a
+    # When residual_threshold=nan there are no inliers and a
     # ValueError with a message should be raised
     base_estimator = LinearRegression()
     ransac_estimator = RANSACRegressor(base_estimator, min_samples=2,
-                                       residual_threshold=0.0, random_state=0,
-                                       max_trials=5)
+                                       residual_threshold=float('nan'),
+                                       random_state=0, max_trials=5)
 
     msg = ("RANSAC could not find a valid consensus set")
     assert_raises_regexp(ValueError, msg, ransac_estimator.fit, X, y)
@@ -538,6 +538,19 @@ def test_ransac_final_model_fit_sample_weight():
     assert_allclose(ransac.estimator_.coef_, final_model.coef_, atol=1e-12)
 
 
+def test_perfect_horizontal_line():
+    # regression test for issue #19497
+    X = np.arange(100)[:, None]
+    y = np.zeros((100, ))
+
+    base_estimator = LinearRegression()
+    ransac_estimator = RANSACRegressor(base_estimator, random_state=0)
+    ransac_estimator.fit(X, y)
+
+    assert_allclose(ransac_estimator.estimator_.coef_, 0.)
+    assert_allclose(ransac_estimator.estimator_.intercept_, 0.)
+
+
 # TODO: Remove in v1.2
 def test_loss_squared_loss_deprecated():
     est1 = RANSACRegressor(loss="squared_loss", random_state=0)


Code review comments:
diff_hunk: @@ -95,7 +95,8 @@ class RANSACRegressor(MetaEstimatorMixin, RegressorMixin,
     residual_threshold : float, default=None
         Maximum residual for a data sample to be classified as an inlier.
         By default the threshold is chosen as the MAD (median absolute
-        deviation) of the target values `y`.
+        deviation) of the target values `y`. Residual strictly equal to the
+        threshold is considered valid.
path: sklearn/linear_model/_ransac.py
text: @user1:
```suggestion
        deviation) of the target values `y`. Points whose residuals are
        strictly equal to the threshold are considered as inliers.
```