Merge pull request scikit-learn#5182 from MechCoder/predict_proba_fix

[MRG + 2] predict_proba should use the softmax function in the multinomial case

Author: GaelVaroquaux

sklearn/linear_model/base.py

@@ -255,7 +255,7 @@ def _predict_proba_lr(self, X):
         np.exp(prob, prob)
         prob += 1
         np.reciprocal(prob, prob)
-        if len(prob.shape) == 1:
+        if prob.ndim == 1:
             return np.vstack([1 - prob, prob]).T
         else:
             # OvR normalization, like LibLinear's predict_probability

sklearn/linear_model/logistic.py

@@ -25,7 +25,7 @@
                              squared_norm)
 from ..utils.optimize import newton_cg
 from ..utils.validation import (as_float_array, DataConversionWarning,
-                                check_X_y)
+                                check_X_y, NotFittedError)
 from ..utils.fixes import expit
 from ..externals.joblib import Parallel, delayed
 from ..cross_validation import check_cv
@@ -1088,6 +1088,13 @@ def predict_proba(self, X):
         The returned estimates for all classes are ordered by the
         label of classes.
 
+        For a multi_class problem, if multi_class is set to be "multinomial"
+        the softmax function is used to find the predicted probability of
+        each class.
+        Else use a one-vs-rest approach, i.e calculate the probability
+        of each class assuming it to be positive using the logistic function.
+        and normalize these values across all the classes.
+
         Parameters
         ----------
         X : array-like, shape = [n_samples, n_features]
@@ -1098,7 +1105,17 @@ def predict_proba(self, X):
             Returns the probability of the sample for each class in the model,
             where classes are ordered as they are in ``self.classes_``.
         """
-        return self._predict_proba_lr(X)
+        if not hasattr(self, "coef_"):
+            raise NotFittedError("Call fit before prediction")
+        calculate_ovr = self.coef_.shape[0] == 1 or self.multi_class == "ovr"
+        if calculate_ovr:
+            return super(LogisticRegression, self)._predict_proba_lr(X)
+        else:
+            prob = self.decision_function(X)
+            np.exp(prob, prob)
+            sum_prob = np.sum(prob, axis=1).reshape((-1, 1))
+            prob /= sum_prob
+            return prob
 
     def predict_log_proba(self, X):
         """Log of probability estimates.

sklearn/linear_model/tests/test_logistic.py

@@ -23,6 +23,7 @@
     )
 from sklearn.cross_validation import StratifiedKFold
 from sklearn.datasets import load_iris, make_classification
+from sklearn.metrics import log_loss
 
 
 X = [[-1, 0], [0, 1], [1, 1]]
@@ -675,3 +676,24 @@ def test_logreg_cv_penalty():
     lr = LogisticRegression(penalty="l1", C=1.0, solver='liblinear')
     lr.fit(X, y)
     assert_equal(np.count_nonzero(lr_cv.coef_), np.count_nonzero(lr.coef_))
+
+
+def test_logreg_predict_proba_multinomial():
+    X, y = make_classification(
+        n_samples=10, n_features=20, random_state=0, n_classes=3, n_informative=10)
+
+    # Predicted probabilites using the true-entropy loss should give a smaller loss
+    # than those using the ovr method.
+    clf_multi = LogisticRegression(multi_class="multinomial", solver="lbfgs")
+    clf_multi.fit(X, y)
+    clf_multi_loss = log_loss(y, clf_multi.predict_proba(X))
+    clf_ovr = LogisticRegression(multi_class="ovr", solver="lbfgs")
+    clf_ovr.fit(X, y)
+    clf_ovr_loss = log_loss(y, clf_ovr.predict_proba(X))
+    assert_greater(clf_ovr_loss, clf_multi_loss)
+
+    # Predicted probabilites using the soft-max function should give a smaller loss
+    # than those using the logistic function.
+    clf_multi_loss = log_loss(y, clf_multi.predict_proba(X))
+    clf_wrong_loss = log_loss(y, clf_multi._predict_proba_lr(X))
+    assert_greater(clf_wrong_loss, clf_multi_loss)