ENH: Add dynamic maximum trial determination to RANSACRegressor

Author: ahojnnes

sklearn/linear_model/ransac.py

@@ -12,6 +12,43 @@
 from .base import LinearRegression
 
 
+_EPSILON = np.spacing(1)
+
+
+def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
+    """Determine number trials such that at least one outlier-free subset is
+    sampled for the given inlier/outlier ratio.
+
+    Parameters
+    ----------
+    n_inliers : int
+        Number of inliers in the data.
+
+    n_samples : int
+        Total number of samples in the data.
+
+    min_samples : int
+        Minimum number of samples chosen randomly from original data.
+
+    probability : float
+        Probability (confidence) that one outlier-free sample is generated.
+
+    Returns
+    -------
+    trials : int
+        Number of trials.
+
+    """
+    inlier_ratio = n_inliers / float(n_samples)
+    nom = max(_EPSILON, 1 - probability)
+    denom = max(_EPSILON, 1 - inlier_ratio ** min_samples)
+    if nom == 1:
+        return 0
+    if denom == 1:
+        return float('inf')
+    return abs(float(np.ceil(np.log(nom) / np.log(denom))))
+
+
 class RANSACRegressor(BaseEstimator, MetaEstimatorMixin, RegressorMixin):
     """RANSAC (RANdom SAmple Consensus) algorithm.
 
@@ -44,7 +81,8 @@ class RANSACRegressor(BaseEstimator, MetaEstimatorMixin, RegressorMixin):
         Minimum number of samples chosen randomly from original data. Treated
         as an absolute number of samples for `min_samples >= 1`, treated as a
         relative number `ceil(min_samples * X.shape[0]`) for
-        `min_samples < 1`. By default a
+        `min_samples < 1`. This is typically chosen as the minimal number of
+        samples necessary to estimate the given `base_estimator`. By default a
         ``sklearn.linear_model.LinearRegression()`` estimator is assumed and
         `min_samples` is chosen as ``X.shape[1] + 1``.
 
@@ -75,6 +113,17 @@ class RANSACRegressor(BaseEstimator, MetaEstimatorMixin, RegressorMixin):
     stop_score : float, optional
         Stop iteration if score is greater equal than this threshold.
 
+    stop_probability : float in range [0, 1], optional
+        RANSAC iteration stops if at least one outlier-free set of the training
+        data is sampled in RANSAC. This requires to generate at least N
+        samples (iterations)::
+
+            N >= log(1 - probability) / log(1 - e**m)
+
+        where the probability (confidence) is typically set to high value such
+        as 0.99 (the default) and e is the current fraction of inliers w.r.t.
+        the total number of samples.
+
     residual_metric : callable, optional
         Metric to reduce the dimensionality of the residuals to 1 for
         multi-dimensional target values ``y.shape[1] > 1``. By default the sum
@@ -110,7 +159,8 @@ def __init__(self, base_estimator=None, min_samples=None,
                  residual_threshold=None, is_data_valid=None,
                  is_model_valid=None, max_trials=100,
                  stop_n_inliers=np.inf, stop_score=np.inf,
-                 residual_metric=None, random_state=None):
+                 stop_probability=0.99, residual_metric=None,
+                 random_state=None):
 
         self.base_estimator = base_estimator
         self.min_samples = min_samples
@@ -120,6 +170,7 @@ def __init__(self, base_estimator=None, min_samples=None,
         self.max_trials = max_trials
         self.stop_n_inliers = stop_n_inliers
         self.stop_score = stop_score
+        self.stop_probability = stop_probability
         self.residual_metric = residual_metric
         self.random_state = random_state
 
@@ -164,6 +215,9 @@ def fit(self, X, y):
             raise ValueError("`min_samples` may not be larger than number "
                              "of samples ``X.shape[0]``.")
 
+        if self.stop_probability < 0 or self.stop_probability > 1:
+            raise ValueError("`stop_probability` must be in range [0, 1].")
+
         if self.residual_threshold is None:
             # MAD (median absolute deviation)
             residual_threshold = np.median(np.abs(y - np.median(y)))
@@ -258,7 +312,11 @@ def fit(self, X, y):
 
             # break if sufficient number of inliers or score is reached
             if (n_inliers_best >= self.stop_n_inliers
-                    or score_best >= self.stop_score):
+                    or score_best >= self.stop_score
+                    or self.n_trials_
+                       >= _dynamic_max_trials(n_inliers_best, n_samples,
+                                              min_samples,
+                                              self.stop_probability)):
                 break
 
         # if none of the iterations met the required criteria

sklearn/linear_model/tests/test_ransac.py

@@ -4,6 +4,7 @@
 
 from sklearn.utils.testing import assert_less
 from sklearn.linear_model import LinearRegression, RANSACRegressor
+from sklearn.linear_model.ransac import _dynamic_max_trials
 
 
 # Generate coordinates of line
@@ -84,7 +85,7 @@ def test_ransac_max_trials():
                                        random_state=0)
     assert getattr(ransac_estimator, 'n_trials_', None) is None
     ransac_estimator.fit(X, y)
-    assert_equal(ransac_estimator.n_trials_, 11)
+    assert_equal(ransac_estimator.n_trials_, 2)
 
 
 def test_ransac_stop_n_inliers():
@@ -277,7 +278,6 @@ def test_ransac_residual_metric():
 
 
 def test_ransac_default_residual_threshold():
-
     base_estimator = LinearRegression()
     ransac_estimator = RANSACRegressor(base_estimator, min_samples=2,
                                        random_state=0)
@@ -293,5 +293,45 @@ def test_ransac_default_residual_threshold():
     assert_equal(ransac_estimator.inlier_mask_, ref_inlier_mask)
 
 
+def test_ransac_dynamic_max_trials():
+    # Numbers hand-calculated and confirmed on page 119 (Table 4.3) in
+    #   Hartley, R.~I. and Zisserman, A., 2004,
+    #   Multiple View Geometry in Computer Vision, Second Edition,
+    #   Cambridge University Press, ISBN: 0521540518
+
+    # e = 0%, min_samples = X
+    assert_equal(_dynamic_max_trials(100, 100, 2, 0.99), 1)
+
+    # e = 5%, min_samples = 2
+    assert_equal(_dynamic_max_trials(95, 100, 2, 0.99), 2)
+    # e = 10%, min_samples = 2
+    assert_equal(_dynamic_max_trials(90, 100, 2, 0.99), 3)
+    # e = 30%, min_samples = 2
+    assert_equal(_dynamic_max_trials(70, 100, 2, 0.99), 7)
+    # e = 50%, min_samples = 2
+    assert_equal(_dynamic_max_trials(50, 100, 2, 0.99), 17)
+
+    # e = 5%, min_samples = 8
+    assert_equal(_dynamic_max_trials(95, 100, 8, 0.99), 5)
+    # e = 10%, min_samples = 8
+    assert_equal(_dynamic_max_trials(90, 100, 8, 0.99), 9)
+    # e = 30%, min_samples = 8
+    assert_equal(_dynamic_max_trials(70, 100, 8, 0.99), 78)
+    # e = 50%, min_samples = 8
+    assert_equal(_dynamic_max_trials(50, 100, 8, 0.99), 1177)
+
+    # e = 0%, min_samples = 10
+    assert_equal(_dynamic_max_trials(1, 100, 10, 0), 0)
+    assert_equal(_dynamic_max_trials(1, 100, 10, 1), float('inf'))
+
+    base_estimator = LinearRegression()
+    ransac_estimator = RANSACRegressor(base_estimator, min_samples=2,
+                                       stop_probability=-0.1)
+    assert_raises(ValueError, ransac_estimator.fit, X, y)
+    ransac_estimator = RANSACRegressor(base_estimator, min_samples=2,
+                                       stop_probability=1.1)
+    assert_raises(ValueError, ransac_estimator.fit, X, y)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()