Repository: scikit-learn/scikit-learn
Pull request: [MRG + 1] Fix the cross_val_predict function for method='predict_proba'
#### Reference Issue
Fixes #7863 

#### What does this implement/fix? Explain your changes.
cross_val_predict() did not specify the order of the classes when method='predict_proba' was passed as the argument. So firstly, the change mentions in the documentation that the columns returned would be sorted. Secondly, for the estimators having a `classes_` attribute, it ensures that the predictions are returned in a sorted manner.

#### Any other comments?
Currently, I haven't yet tested the code. Would be completing that within a day or so.



Problem statement:
order of returned probabilites unclear for cross_val_predict with method=predict_proba
Cross_val_predict has a new method parameter which is typically set to "predict_proba" to return probabilities for each class.

However the order of the classes returned is unclear. Either self.classes_ needs to be set; or the results need to be returned in a predictable order. Otherwise we have a list of probabilities for each class but no way to know which column relates to which class.

Hints from issue discussion:
Yes, I suppose this (and disagreements in the set of classes between splits) was overlooked.

(For internal models, I'm fairly sure that `classes_` is always alphabetically ordered.)

Can I start working on this?

After running a few tests on the iris dataset using LogisticRegression as the estimator, it became clear that the order of the classes appearing did not matter in the final result. [This](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression.predict_log_proba) gives a simple and clear explanation. 

Iris aside, it's possible to create a cross-validation strategy that will
land up with only a subset of the classes being present in the training
data for a particular training set. This should be handled.

Firstly the documentation should point out that the classes will be in
sorted order, and this should, in theory at least, be confirmed, reordering
if necessary, from the underlying estimators.

On 16 November 2016 at 11:41, Aman Dalmia notifications@github.com wrote:

> After running a few tests on the iris datasets using LogisticRegression as
> the estimator, it became clear that the order of the classes appearing did
> not matter in the final result. This
> http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html#sklearn.linear_model.LogisticRegression.predict_log_proba
> gives a simple and clear explanation.
> 
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> https://github.com/scikit-learn/scikit-learn/issues/7863#issuecomment-260817612,
> or mute the thread
> https://github.com/notifications/unsubscribe-auth/AAEz65Z6-gNAXqLB0OC8OTlEeBKG9pBeks5q-lFNgaJpZM4Kw2SU
> .

True, we could miss out on a few observations in cross-validation. As for the confirmation of the classes being sorted, since these lines of the cross_val_predict function already ensure that only estimators capable of calculating the probability are passed 'predict_proba' as the method, and all such estimators return the classes in sorted order, do we need any other mode of confirmation?
Please correct me if I misunderstood what you intended to say.

``` python
 if not callable(getattr(estimator, method)):
        raise AttributeError('{} not implemented in estimator'
                             .format(method))
```

Well, we don't promise that all return classes in sorted order; all store
classes in `.classes_`

On 16 November 2016 at 14:44, Aman Dalmia notifications@github.com wrote:

> True, we could miss out on a few observations in cross-validation. As for
> the confirmation of the classes being sorted, since these lines of the
> cross_val_predict function already ensure that only estimators capable of
> calculating the probability are passed 'predict_proba' as the method, and
> all such estimators return the classes in sorted order, do we need any
> other mode of confirmation?
> Please correct me if I misunderstood what you intended to say.
> 
>  if not callable(getattr(estimator, method)):
>         raise AttributeError('{} not implemented in estimator'
>                              .format(method))
> 
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> https://github.com/scikit-learn/scikit-learn/issues/7863#issuecomment-260845838,
> or mute the thread
> https://github.com/notifications/unsubscribe-auth/AAEz643iTLmmU5CBwLaMK5h6Iwqak0f3ks5q-nwFgaJpZM4Kw2SU
> .

Then how about, in the `_fit_and_predict` function, I retrieve the `classes_` attribute, make a dictionary of the column indices and the class labels, sort them on the class labels and apply the transformation to the predictions, thereby ensuring that a sorted result is returned?

That sort of thing, yes

On 16 November 2016 at 16:23, Aman Dalmia notifications@github.com wrote:

> Then how about, in the _fit_and_predict function, I retrieve the 'classes_'
> attribute, make a dictionary of the column indices and the class labels,
> sort them and apply the transformation to the predictions, thereby ensuring
> that a sorted result is returned?
> 
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> https://github.com/scikit-learn/scikit-learn/issues/7863#issuecomment-260857181,
> or mute the thread
> https://github.com/notifications/unsubscribe-auth/AAEz61TyLjZQJX7Jjehj00p84uYkIiMHks5q-pNWgaJpZM4Kw2SU
> .

That would be a good idea as it would make sure 3rd party classifiers used the same standard. I know it works today for LogisticRegression and RandomForestClassifier. However I am unsure for XGBClassifier as that is outside sklearn.

Commit message: FIX: fixed _fit_and_predict and added tests

Patch to review:
diff --git a/sklearn/model_selection/_validation.py b/sklearn/model_selection/_validation.py
index 23db2a9cebc77..329f7308988a3 100644
--- a/sklearn/model_selection/_validation.py
+++ b/sklearn/model_selection/_validation.py
@@ -365,7 +365,9 @@ def cross_val_predict(estimator, X, y=None, groups=None, cv=None, n_jobs=1,
               as in '2*n_jobs'
 
     method : string, optional, default: 'predict'
-        Invokes the passed method name of the passed estimator.
+        Invokes the passed method name of the passed estimator. For
+        method='predict_proba', the columns correspond to the classes
+        in sorted order.
 
     Returns
     -------
@@ -474,6 +476,18 @@ def _fit_and_predict(estimator, X, y, train, test, verbose, fit_params,
         estimator.fit(X_train, y_train, **fit_params)
     func = getattr(estimator, method)
     predictions = func(X_test)
+    if method in ['decision_function', 'predict_proba', 'predict_log_proba']:
+        true_classes = np.unique(y)
+        train_classes = np.unique(y_train)
+        predictions_ = np.zeros((X_test.shape[0], true_classes.shape[0]))
+        if method is 'decision_function' and len(train_classes) == 2:
+            class_predictions = estimator.predict(X_test)
+            for i, j in enumerate(class_predictions):
+                predictions_[i, j] = predictions[i]
+        else:
+            for i, j in enumerate(train_classes):
+                predictions_[:, j] = predictions[:, i]
+        predictions = predictions_
     return predictions, test
 
 
diff --git a/sklearn/model_selection/tests/test_validation.py b/sklearn/model_selection/tests/test_validation.py
index 31c5fc8257528..6b1b81350f772 100644
--- a/sklearn/model_selection/tests/test_validation.py
+++ b/sklearn/model_selection/tests/test_validation.py
@@ -914,6 +914,10 @@ def test_cross_val_predict_with_method():
     X, y = shuffle(X, y, random_state=0)
     classes = len(set(y))
 
+    y[:50] = 0
+    y[50:100] = 1
+    X[:100], y[:100] = shuffle(X[:100], y[:100], random_state=0)
+
     kfold = KFold(len(iris.target))
 
     methods = ['decision_function', 'predict_proba', 'predict_log_proba']
@@ -929,7 +933,15 @@ def test_cross_val_predict_with_method():
         # Naive loop (should be same as cross_val_predict):
         for train, test in kfold.split(X, y):
             est.fit(X[train], y[train])
-            expected_predictions[test] = func(X[test])
+            train_classes = np.unique(y[train])
+            expected_predictions_ = func(X[test])
+            if method is 'decision_function' and len(train_classes) == 2:
+                expected_test_predictions = est.predict(X[test])
+                for i, j in enumerate(expected_test_predictions):
+                    expected_predictions[i, j] = expected_predictions_[i]
+            else:
+                for i, j in enumerate(train_classes):
+                    expected_predictions[test, j] = expected_predictions_[:, i]
 
         predictions = cross_val_predict(est, X, y, method=method,
                                         cv=kfold)


Code review comments:
diff_hunk: @@ -474,6 +476,18 @@ def _fit_and_predict(estimator, X, y, train, test, verbose, fit_params,
         estimator.fit(X_train, y_train, **fit_params)
     func = getattr(estimator, method)
     predictions = func(X_test)
+    if method in ['decision_function', 'predict_proba', 'predict_log_proba']:
+        true_classes = np.unique(y)
+        train_classes = np.unique(y_train)
path: sklearn/model_selection/_validation.py
text: @user1:
I think we should be using the estimator's `classes_` here

@author:
But a few estimators implementing the `predict_proba` method, don't seem to have `classes_` as an attribute like `SVC`.

@user1:
Huh?
```python
>>> SVC().fit(np.arange(20).reshape(5, 4), [0, 1, 1, 0, 1]).classes_
array([0, 1])
```

@author:
Sorry, I just looked at the [documentation](http://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html) where it's not mentioned in the `Attributes` section and made the note. I'll add the change.

diff_hunk: @@ -914,6 +914,10 @@ def test_cross_val_predict_with_method():
     X, y = shuffle(X, y, random_state=0)
     classes = len(set(y))
 
+    y[:50] = 0
path: sklearn/model_selection/tests/test_validation.py
text: @user1:
It would be more normal to make a separate test for tricky cases, rather than overload the test for basic functionality. This helps diagnose, when a test fails, what went wrong.

Certainly, this dataset modification needs a comment as to its purpose.

diff_hunk: @@ -474,6 +476,18 @@ def _fit_and_predict(estimator, X, y, train, test, verbose, fit_params,
         estimator.fit(X_train, y_train, **fit_params)
     func = getattr(estimator, method)
     predictions = func(X_test)
+    if method in ['decision_function', 'predict_proba', 'predict_log_proba']:
+        true_classes = np.unique(y)
+        train_classes = np.unique(y_train)
+        predictions_ = np.zeros((X_test.shape[0], true_classes.shape[0]))
+        if method is 'decision_function' and len(train_classes) == 2:
+            class_predictions = estimator.predict(X_test)
path: sklearn/model_selection/_validation.py
text: @user1:
I don't see what the output of `predict` has to do with it?

@author:
I realize I made a mistake here. Will correct this.

diff_hunk: @@ -474,6 +476,18 @@ def _fit_and_predict(estimator, X, y, train, test, verbose, fit_params,
         estimator.fit(X_train, y_train, **fit_params)
     func = getattr(estimator, method)
     predictions = func(X_test)
+    if method in ['decision_function', 'predict_proba', 'predict_log_proba']:
+        true_classes = np.unique(y)
+        train_classes = np.unique(y_train)
+        predictions_ = np.zeros((X_test.shape[0], true_classes.shape[0]))
+        if method is 'decision_function' and len(train_classes) == 2:
+            class_predictions = estimator.predict(X_test)
+            for i, j in enumerate(class_predictions):
+                predictions_[i, j] = predictions[i]
+        else:
+            for i, j in enumerate(train_classes):
path: sklearn/model_selection/_validation.py
text: @user1:
I don't see how this is mapping into the `true_classes` space

@author:
With the help of the example that I used in my test case, `y` has the classes`[0, 1, 2]`. If one block of `kfold.split`, `y_train` only has classes `[0,2]`, then this replaces the corresponding columns in `_predictions` (50 x 3) with the appropriate columns from `predictions` (50 x 2) and leaves the other columns as zeros.

diff_hunk: @@ -474,6 +476,18 @@ def _fit_and_predict(estimator, X, y, train, test, verbose, fit_params,
         estimator.fit(X_train, y_train, **fit_params)
     func = getattr(estimator, method)
     predictions = func(X_test)
+    if method in ['decision_function', 'predict_proba', 'predict_log_proba']:
path: sklearn/model_selection/_validation.py
text: @user1:
I wonder if we should be inspecting the shape of the prediction, rather than doing this based on name.

@author:
The problem lies is that `predict_proba` and `predict_log_proba` return a 2D array when there are 2 classes in the training set but `decision_function` returns a 1D array. So, I chose to go by method to differentiate the two.