Repository: scikit-learn/scikit-learn
Pull request: fixed a bug in isotonic regression
#### Reference Issue
Fixes #9432 

#### What does this implement/fix? Explain your changes.
There was an apparent bug in the preprocessing steps for isotonic regression (_make_unique), in that the weights for duplicate points were averaged, rather that summed, effectively down-weighting duplicate points in the actual isotonic regression. This fix changes the default behaviour to summing the weights, such that isotonic regression now returns the correct predictions. However, I also implemented an option (do_weight_averaging; default=False) to preserve the previous behaviour as an option.

#### Any other comments?
Here is a simple example to demonstrate the change:

```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.isotonic import IsotonicRegression
x = np.array([0, 1, 1, 2, 3])
y = np.array([0, 0, 1, 0, 1])
weights = np.ones(len(y))
ir = IsotonicRegression()
ir.fit(x, y, sample_weight=weights)
y_pred = ir.predict(x)
print("Predictions: %s" % str(y_pred))
squared_error = np.dot(weights, (y - y_pred)**2)
print("squared error = %0.3f" % squared_error)
```

Under the old code, this produced:
> Predictions: [ 0.    0.25  0.25  0.25  1.  ]
> squared error = 0.688

After the update, this produces:
> Predictions: [ 0.          0.33333333  0.33333333  0.33333333  1.        ]
> squared error = 0.667

Setting do_weight_averaging to true (i.e. ir = IsotonicRegression(do_weight_averaging=True)) restores the original result.


Problem statement:
Incorrect predictions from isotonic regression in the presence of repeated X values
#### Description
Isotonic regression returns an incorrect solution in the presence of repeated X values. 

The problem seems to be that fit() begins by taking the weighted average of points with the same x value. However, rather than associating the resulting average point with the sum of the weights, it is using the average of the weights (in _make_unique(), in  _isotonic.pyx, lines 103 and 116). When the fitting the isotonic regression, the combined points will then count for less than they should (according to the objective), and incorrect fitted values are returned.

It is conceivable that this is the intended behaviour (i.e. if the objective in the documentation refers to the output of _make_unique, rather than the original data), but it seems to conflict with the documentation, which explicitly refers to the possibility of equality among X values. The docs also mention that ties are broken using the "secondary method from Leeuw, 1977", but that too seems to sum the weights rather than averaging them (page 2, 6th line from the bottom). 

If I have missed something, and this is the intended behaviour, I think we should clarify the documentation, and possibly add an option to not average in this way. However, I would advocate for summing rather than averaging as the default.

#### Steps/Code to Reproduce
This code presents a simple example that returns an incorrect result (according to the stated objective). It then uses the sample_weight option to artificially up-weight the points which get combined (at x=1). This, results in a lower squared error on the original data, demonstrating that the default behaviour does not achieve the minimum. 

```
import numpy as np
import matplotlib.pyplot as plt
from sklearn.isotonic import IsotonicRegression

x = np.array([0, 1, 1, 2, 3])
y = np.array([0, 0, 1, 0, 1])

ir = IsotonicRegression()
ir.fit(x, y)
y_pred = ir.predict(x)

squared_error = np.sum((y - y_pred)**2)
print("Default SE = %0.3f" % squared_error)

fig, ax = plt.subplots()
ax.scatter(x, y)
ax.plot(x, y_pred, label="default")

# force IR to do the right thing by artificially up-weighting the data
sample_weight = np.array([1.0, 2.0, 2.0, 1.0, 1.0])
ir.fit(x, y, sample_weight=sample_weight)
y_pred = ir.predict(x)

squared_error = np.sum((y - y_pred)**2)
print("Reweighted SE = %0.3f" % squared_error)

ax.plot(x, y_pred, label="reweighted")
ax.legend()
plt.show()
```

#### Expected Results
I believe the output of isotonic regression should predict [0, 0.33, 0.33, 0.33, 1] at the original input points, with squared error 0.667.  _make_unique should produce x = (0, 1, 2, 3), y = (0, 0.5, 0, 1), sample_weight = (1, 2, 1, 1).

#### Actual Results
The actual predicted values at these points are [0, 0.25, 0.25, 0.25, 1], with squared error 0.688. _make_unique is producing x = (0, 1, 2, 3), y = (0, 0.5, 0, 1), sample_weight = (1, 1, 1, 1).

#### Versions
Darwin-16.6.0-x86_64-i386-64bit
Python 3.6.1 |Continuum Analytics, Inc.| (default, Mar 22 2017, 19:25:17) 
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]
NumPy 1.12.1
SciPy 0.19.0
Scikit-Learn 0.18.2

Hints from issue discussion:

Commit message: added non-regression test

Patch to review:
diff --git a/sklearn/_isotonic.pyx b/sklearn/_isotonic.pyx
index 1cec075fc6fc7..55f3bea7ac6d5 100644
--- a/sklearn/_isotonic.pyx
+++ b/sklearn/_isotonic.pyx
@@ -70,7 +70,8 @@ def _inplace_contiguous_isotonic_regression(DOUBLE[::1] y, DOUBLE[::1] w):
 @cython.cdivision(True)
 def _make_unique(np.ndarray[dtype=np.float64_t] X,
                  np.ndarray[dtype=np.float64_t] y,
-                 np.ndarray[dtype=np.float64_t] sample_weights):
+                 np.ndarray[dtype=np.float64_t] sample_weights,
+                 do_weight_averaging):
     """Average targets for duplicate X, drop duplicates.
 
     Aggregates duplicate X values into a single X value where
@@ -78,6 +79,10 @@ def _make_unique(np.ndarray[dtype=np.float64_t] X,
     targets.
 
     Assumes that X is ordered, so that all duplicates follow each other.
+
+    If do_weight_averaging is True, the weights for duplicate items will
+    be averaged, effectively down-weighting the duplicates points. If
+    False, the weights for duplicate points will be summed.
     """
     unique_values = len(np.unique(X))
     if unique_values == len(X):
@@ -95,12 +100,16 @@ def _make_unique(np.ndarray[dtype=np.float64_t] X,
     cdef int j
     cdef np.float64_t x
     cdef int n_samples = len(X)
+    cdef bint weight_averaging = do_weight_averaging
     for j in range(n_samples):
         x = X[j]
         if x != current_x:
             # next unique value
             x_out[i] = current_x
-            weights_out[i] = current_weight / current_count
+            if weight_averaging:
+                weights_out[i] = current_weight / current_count
+            else:
+                weights_out[i] = current_weight
             y_out[i] = current_y / current_weight
             i += 1
             current_x = x
@@ -113,6 +122,9 @@ def _make_unique(np.ndarray[dtype=np.float64_t] X,
             current_count += 1
 
     x_out[i] = current_x
-    weights_out[i] = current_weight / current_count
+    if weight_averaging:
+        weights_out[i] = current_weight / current_count
+    else:
+        weights_out[i] = current_weight
     y_out[i] = current_y / current_weight
     return x_out, y_out, weights_out
diff --git a/sklearn/isotonic.py b/sklearn/isotonic.py
index 245fc95ff69b2..be930c71e8d1e 100644
--- a/sklearn/isotonic.py
+++ b/sklearn/isotonic.py
@@ -179,6 +179,11 @@ class IsotonicRegression(BaseEstimator, TransformerMixin, RegressorMixin):
         set to the value corresponding to the nearest train interval endpoint.
         When set to "raise", allow ``interp1d`` to throw ValueError.
 
+    do_weight_averaging : boolean, optional, default: False
+        The ``do_weight_averaging'' parameter determines how weights for points
+        with duplicate X values are handled. If True, the weights for
+        duplicate points will be averaged before fitting the isotonic
+        regression. If False, the weights will be summed.
 
     Attributes
     ----------
@@ -211,11 +216,12 @@ class IsotonicRegression(BaseEstimator, TransformerMixin, RegressorMixin):
     Leeuw, Psychometrica, 1977
     """
     def __init__(self, y_min=None, y_max=None, increasing=True,
-                 out_of_bounds='nan'):
+                 out_of_bounds='nan', do_weight_averaging=False):
         self.y_min = y_min
         self.y_max = y_max
         self.increasing = increasing
         self.out_of_bounds = out_of_bounds
+        self.do_weight_averaging = do_weight_averaging
 
     @property
     @deprecated("Attribute ``X_`` is deprecated in version 0.18 and will be"
@@ -293,7 +299,7 @@ def _build_y(self, X, y, sample_weight, trim_duplicates=True):
         X, y, sample_weight = [array[order].astype(np.float64, copy=False)
                                for array in [X, y, sample_weight]]
         unique_X, unique_y, unique_sample_weight = _make_unique(
-            X, y, sample_weight)
+            X, y, sample_weight, self.do_weight_averaging)
 
         # Store _X_ and _y_ to maintain backward compat during the deprecation
         # period of X_ and y_
diff --git a/sklearn/tests/test_isotonic.py b/sklearn/tests/test_isotonic.py
index d5d0715a0fb7f..ccfb45f50862c 100644
--- a/sklearn/tests/test_isotonic.py
+++ b/sklearn/tests/test_isotonic.py
@@ -166,6 +166,30 @@ def test_isotonic_regression_ties_secondary_():
     assert_array_almost_equal(ir.fit_transform(x, y), y_true, 4)
 
 
+def test_isotonic_regression_with_ties_in_differently_sized_groups():
+    """
+    Non-regression test to handle issue 9432:
+    https://github.com/scikit-learn/scikit-learn/issues/9432
+
+    Compare against output in R:
+    > library("isotone")
+    > x <- c(0, 1, 1, 2, 3, 4)
+    > y <- c(0, 0, 1, 0, 0, 1)
+    > res1 <- gpava(x, y, ties="secondary")
+    > res1$x
+
+    `isotone` version: 1.1-0, 2015-07-24
+    R version: R version 3.3.2 (2016-10-31)
+    """
+    x = np.array([0, 1, 1, 2, 3, 4])
+    y = np.array([0, 0, 1, 0, 0, 1])
+    y_true = np.array([0., 0.25, 0.25, 0.25, 0.25, 1.])
+    ir = IsotonicRegression()
+    ir.fit(x, y)
+    assert_array_almost_equal(ir.transform(x), y_true, 4)
+    assert_array_almost_equal(ir.fit_transform(x, y), y_true, 4)
+
+
 def test_isotonic_regression_reversed():
     y = np.array([10, 9, 10, 7, 6, 6.1, 5])
     y_ = IsotonicRegression(increasing=False).fit_transform(


Code review comments:
diff_hunk: @@ -70,14 +70,19 @@ def _inplace_contiguous_isotonic_regression(DOUBLE[::1] y, DOUBLE[::1] w):
 @cython.cdivision(True)
 def _make_unique(np.ndarray[dtype=np.float64_t] X,
                  np.ndarray[dtype=np.float64_t] y,
-                 np.ndarray[dtype=np.float64_t] sample_weights):
+                 np.ndarray[dtype=np.float64_t] sample_weights,
+                 do_weight_averaging):
path: sklearn/_isotonic.pyx
text: @user1:
Can't you make this a bint directly?