Repository: scikit-learn/scikit-learn
Pull request:  [MRG] FIX solve consistency between predict and predict_proba in AdaBoost 
closes #14084 
closes #2974

compute the probabilities in AdaBoostClassifier as specified in "Multiclass AdaBoost"

Problem statement:
AdaBoost's "SAMME" algorithm uses 'predict' while fitting and 'predict_proba' while predicting probas
Subj. This seems to me to be a wrong approach, moreover this drives to such mistakes:

<pre>
AdaBoostClassifier(algorithm="SAMME", base_estimator=SVC()).fit(trainX, trainY).predict_proba(testX)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-108-1d666912dada> in <module>()
----> 1 AdaBoostClassifier(algorithm="SAMME", base_estimator=SVC()).fit(trainX, trainY).predict_proba(testX)

/Library/Python/2.7/site-packages/sklearn/ensemble/weight_boosting.pyc in predict_proba(self, X)
    716             proba = sum(estimator.predict_proba(X) * w
    717                         for estimator, w in zip(self.estimators_,
--> 718                                                 self.estimator_weights_))
    719 
    720         proba /= self.estimator_weights_.sum()

/Library/Python/2.7/site-packages/sklearn/ensemble/weight_boosting.pyc in <genexpr>((estimator, w))
    715         else:   # self.algorithm == "SAMME"
    716             proba = sum(estimator.predict_proba(X) * w
--> 717                         for estimator, w in zip(self.estimators_,
    718                                                 self.estimator_weights_))
    719 

/Library/Python/2.7/site-packages/sklearn/svm/base.pyc in predict_proba(self, X)
    493         if not self.probability:
    494             raise NotImplementedError(
--> 495                 "probability estimates must be enabled to use this method")
    496 
    497         if self._impl not in ('c_svc', 'nu_svc'):

NotImplementedError: probability estimates must be enabled to use this method
</pre>


AdaboostClassfier predict does not always correspond to class with highest probability
#### Description
AdaboostClassifier `model.predict()` is not always equal to `np.argmax(model.predict_proba(), axis=1)`.

#### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
from sklearn.ensemble import AdaBoostClassifier

data = load_digits()
X = data.data
y = data.target

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=42)
model = AdaBoostClassifier(base_estimator=None, n_estimators=20, learning_rate=1,
                           algorithm='SAMME', random_state=42).fit(X_train, y_train)


print(np.mean(np.argmax(model.predict_proba(X_test), axis=1) == model.predict(X_test)))
```

#### Expected Results
1.0

```python
print(model.predict_proba(X_test)[0]) gives the class probability scores for the first test example:
[0.09949952 0.1000983  0.10002473 0.09970723 0.10019514 0.10022755
 **0.10087143** 0.09974335 0.10005309 0.09957965]
```

Clearly class 6 has the highest probability of 0.10087143, but predict() gives class 5 as the result.

```python
print(model.predict(X_test)[0])
Result: 5
```

#### Actual Results
0.44048943270300334

#### Versions
System:
    python: 3.6.8 (v3.6.8:3c6b436a57, Dec 24 2018, 02:04:31)  [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]
executable: /Users/prroy/Documents/MachineLearning/onnx_projects/onnx_env/bin/python3.6
   machine: Darwin-18.6.0-x86_64-i386-64bit

BLAS:
    macros: NO_ATLAS_INFO=3, HAVE_CBLAS=None
  lib_dirs: 
cblas_libs: cblas

Python deps:
       pip: 19.1.1
setuptools: 41.0.1
   sklearn: 0.21.2
     numpy: 1.16.4
     scipy: 1.1.0
    Cython: None
    pandas: 0.23.4

Hints from issue discussion:
(Not an AdaBoost expert)

Why is it wrong? How else would you define `predict_proba`?

The idea of using only predictions during training and use afterwards probas of base_estimators is strange. The base_estimator can return -0.1 and 0.9 or -0.9 and 0.1.

They will have same predictions and different probas - but you don't take it into account.

The 'standart' scheme as I understand is:
There is score:

score(obj)  = sum[ weight_i \* esimator_i.predict(obj) ],
assuming that predict returns 1 and -1

Then this score is turned to proba by some sigmoid function (score_to_proba in sklearn's gradientBoosting)

> The base_estimator can return -0.1 and 0.9 or -0.9 and 0.1.

Not from `predict`, that returns discrete labels. But it is a bit strange that `predict_proba` should be needed on the base estimator if it's not used in training... @ndawe?

For sure not from predict, sorry. I wanted to say, the predict_proba 
can be [0.1, 0.9] or [0.6, 0.4] of one estimator and
[0.9, 0.1] or [0.4, 0.6] of another, but their predicts will be similar.

This estimators will be considered as similar during training, but at this moment they will have totally different influence on result of predict_proba of AdaBoost

In fact, I don't think `predict_proba` should even be defined when AdaBoost is built from the SAMME variant (which is a discrete boosting algorithm, for base estimators that only support crisp predictions). If I remember correctly, we added that for convenience only. 

What do you think @ndawe?

SAMME.R uses the class probabilities in the training but SAMME does not. That is how those algorithms are designed.

I agree that in the SAMME case, `predict_proba` could be a little ambiguous. Yes, you could transform the discrete labels into some form of probability as you suggest, but the current implementation uses the underlying `predict_proba` of the base estimator. I don't think this is strange but I'm open to suggestions. If the base estimator supports class probabilities, then SAMME.R is the better algorithm (generally). I suppose what we want here is some way of extracting "probabilities" from a boosted model that can only deliver discrete labels.

> I suppose what we want here is some way of extracting "probabilities" from a boosted model that can only deliver discrete labels.

Right. I was working on implementation of uBoost, some variation of AdaBoost, and this was the issue - the predictions of probabilities on some stage are used there to define weights on next iterations.

Somehow that resulted in poor uniformity of predictions. Fixing the `predict_proba` resolved this issue

The change I propose in predict_proba is (if `predict` of estimator returns only 0, 1!)

<pre>
 score = sum((2*estimator.predict(X) - 1) * w
       for estimator, w in zip(self.estimators_,  self.estimator_weights_))
 proba = sigmoid(score)
</pre>

where sigmoid is some sigmoid function.

Looks like some numeric issues, thanks for reporting.
I gave a look at it. It could also be related to #2974
From my understanding, I think that the way we compute the probability (at least with SAMME) is not as I understand it from the paper: https://web.stanford.edu/~hastie/Papers/samme.pdf

Following the eq. (4) of the paper, it seems that we can compute the probabilities from the decision function:

```python
from sklearn.utils.extmath import softmax
y_pred = model.decision_function(X_test)
y_pred *= 1 / (len(model.classes_) - 1)
y_pred = softmax(y_pred, copy=False)
```

This would give us:

```python
In [1]: np.mean(np.argmax(y_pred, axis=1) == model.predict(X_test))
Out [1]: 1.0
In [2]: y_pred.sum(axis=1)
Out [2]: [1., 1., ... 1., 1.]
```

Commit message: DOC add whats new

Patch to review:
diff --git a/doc/whats_new/v0.22.rst b/doc/whats_new/v0.22.rst
index 85b8c330bd9fe..04fdf6a510e66 100644
--- a/doc/whats_new/v0.22.rst
+++ b/doc/whats_new/v0.22.rst
@@ -55,6 +55,11 @@ Changelog
   validation data separately to avoid any data leak. :pr:`13933` by
   `Nicolas Hug`_.
 
+- |Fix| :class:`ensemble.AdaBoostClassifier` computes probabilities based on
+  the decision function as in the literature. Thus, `predict` and
+  `predict_proba` give consistent results.
+  :pr:`14114` by :user:`Guillaume Lemaitre <glemaitre>`.
+
 :mod:`sklearn.linear_model`
 ...........................
 
diff --git a/sklearn/ensemble/tests/test_weight_boosting.py b/sklearn/ensemble/tests/test_weight_boosting.py
index 5b1c38e16d00e..d9e845ce4782c 100755
--- a/sklearn/ensemble/tests/test_weight_boosting.py
+++ b/sklearn/ensemble/tests/test_weight_boosting.py
@@ -1,6 +1,7 @@
 """Testing for the boost module (sklearn.ensemble.boost)."""
 
 import numpy as np
+import pytest
 
 from sklearn.utils.testing import assert_array_equal, assert_array_less
 from sklearn.utils.testing import assert_array_almost_equal
@@ -84,15 +85,15 @@ def test_oneclass_adaboost_proba():
     assert_array_almost_equal(clf.predict_proba(X), np.ones((len(X), 1)))
 
 
-def test_classification_toy():
+@pytest.mark.parametrize("algorithm", ["SAMME", "SAMME.R"])
+def test_classification_toy(algorithm):
     # Check classification on a toy dataset.
-    for alg in ['SAMME', 'SAMME.R']:
-        clf = AdaBoostClassifier(algorithm=alg, random_state=0)
-        clf.fit(X, y_class)
-        assert_array_equal(clf.predict(T), y_t_class)
-        assert_array_equal(np.unique(np.asarray(y_t_class)), clf.classes_)
-        assert_equal(clf.predict_proba(T).shape, (len(T), 2))
-        assert_equal(clf.decision_function(T).shape, (len(T),))
+    clf = AdaBoostClassifier(algorithm=algorithm, random_state=0)
+    clf.fit(X, y_class)
+    assert_array_equal(clf.predict(T), y_t_class)
+    assert_array_equal(np.unique(np.asarray(y_t_class)), clf.classes_)
+    assert clf.predict_proba(T).shape == (len(T), 2)
+    assert clf.decision_function(T).shape == (len(T),)
 
 
 def test_regression_toy():
@@ -151,32 +152,31 @@ def test_boston():
                  len(reg.estimators_))
 
 
-def test_staged_predict():
+@pytest.mark.parametrize("algorithm", ["SAMME", "SAMME.R"])
+def test_staged_predict(algorithm):
     # Check staged predictions.
     rng = np.random.RandomState(0)
     iris_weights = rng.randint(10, size=iris.target.shape)
     boston_weights = rng.randint(10, size=boston.target.shape)
 
-    # AdaBoost classification
-    for alg in ['SAMME', 'SAMME.R']:
-        clf = AdaBoostClassifier(algorithm=alg, n_estimators=10)
-        clf.fit(iris.data, iris.target, sample_weight=iris_weights)
+    clf = AdaBoostClassifier(algorithm=algorithm, n_estimators=10)
+    clf.fit(iris.data, iris.target, sample_weight=iris_weights)
 
-        predictions = clf.predict(iris.data)
-        staged_predictions = [p for p in clf.staged_predict(iris.data)]
-        proba = clf.predict_proba(iris.data)
-        staged_probas = [p for p in clf.staged_predict_proba(iris.data)]
-        score = clf.score(iris.data, iris.target, sample_weight=iris_weights)
-        staged_scores = [
-            s for s in clf.staged_score(
-                iris.data, iris.target, sample_weight=iris_weights)]
-
-        assert_equal(len(staged_predictions), 10)
-        assert_array_almost_equal(predictions, staged_predictions[-1])
-        assert_equal(len(staged_probas), 10)
-        assert_array_almost_equal(proba, staged_probas[-1])
-        assert_equal(len(staged_scores), 10)
-        assert_array_almost_equal(score, staged_scores[-1])
+    predictions = clf.predict(iris.data)
+    staged_predictions = [p for p in clf.staged_predict(iris.data)]
+    proba = clf.predict_proba(iris.data)
+    staged_probas = [p for p in clf.staged_predict_proba(iris.data)]
+    score = clf.score(iris.data, iris.target, sample_weight=iris_weights)
+    staged_scores = [
+        s for s in clf.staged_score(
+            iris.data, iris.target, sample_weight=iris_weights)]
+
+    assert_equal(len(staged_predictions), 10)
+    assert_array_almost_equal(predictions, staged_predictions[-1])
+    assert_equal(len(staged_probas), 10)
+    assert_array_almost_equal(proba, staged_probas[-1])
+    assert_equal(len(staged_scores), 10)
+    assert_array_almost_equal(score, staged_scores[-1])
 
     # AdaBoost regression
     clf = AdaBoostRegressor(n_estimators=10, random_state=0)
@@ -505,3 +505,20 @@ def test_multidimensional_X():
     boost = AdaBoostRegressor(DummyRegressor())
     boost.fit(X, yr)
     boost.predict(X)
+
+
+@pytest.mark.parametrize("algorithm", ["SAMME", "SAMME.R"])
+def test_adaboost_consistent_predict(algorithm):
+    # check that predict_proba and predict give consistent results
+    # regression test for:
+    # https://github.com/scikit-learn/scikit-learn/issues/14084
+    X_train, X_test, y_train, y_test = train_test_split(
+        *datasets.load_digits(return_X_y=True), random_state=42
+    )
+    model = AdaBoostClassifier(algorithm=algorithm, random_state=42)
+    model.fit(X_train, y_train)
+
+    assert_array_equal(
+        np.argmax(model.predict_proba(X_test), axis=1),
+        model.predict(X_test)
+    )
diff --git a/sklearn/ensemble/weight_boosting.py b/sklearn/ensemble/weight_boosting.py
index 3cb4baa0d9a0c..3b3c7b819c2d2 100644
--- a/sklearn/ensemble/weight_boosting.py
+++ b/sklearn/ensemble/weight_boosting.py
@@ -34,6 +34,7 @@
 
 from ..tree import DecisionTreeClassifier, DecisionTreeRegressor
 from ..utils import check_array, check_random_state, check_X_y, safe_indexing
+from ..utils.extmath import softmax
 from ..utils.extmath import stable_cumsum
 from ..metrics import accuracy_score, r2_score
 from ..utils.validation import check_is_fitted
@@ -748,6 +749,24 @@ class in ``classes_``, respectively.
             else:
                 yield pred / norm
 
+    @staticmethod
+    def _compute_proba_from_decision(decision, n_classes):
+        """Compute probabilities from the decision function.
+
+        This is based eq. (4) of [1] where:
+            p(y=c|X) = exp((1 / K-1) f_c(X)) / sum_k(exp((1 / K-1) f_k(X)))
+                     = softmax((1 / K-1) * f(X))
+
+        References
+        ----------
+        .. [1] J. Zhu, H. Zou, S. Rosset, T. Hastie, "Multi-class AdaBoost",
+               2009.
+        """
+        if n_classes == 2:
+            decision = np.vstack([-decision, decision]).T
+        decision *= 1 / (n_classes - 1)
+        return softmax(decision, copy=False)
+
     def predict_proba(self, X):
         """Predict class probabilities for X.
 
@@ -775,22 +794,8 @@ def predict_proba(self, X):
         if n_classes == 1:
             return np.ones((_num_samples(X), 1))
 
-        if self.algorithm == 'SAMME.R':
-            # The weights are all 1. for SAMME.R
-            proba = sum(_samme_proba(estimator, n_classes, X)
-                        for estimator in self.estimators_)
-        else:  # self.algorithm == "SAMME"
-            proba = sum(estimator.predict_proba(X) * w
-                        for estimator, w in zip(self.estimators_,
-                                                self.estimator_weights_))
-
-        proba /= self.estimator_weights_.sum()
-        proba = np.exp((1. / (n_classes - 1)) * proba)
-        normalizer = proba.sum(axis=1)[:, np.newaxis]
-        normalizer[normalizer == 0.0] = 1.0
-        proba /= normalizer
-
-        return proba
+        decision = self.decision_function(X)
+        return self._compute_proba_from_decision(decision, n_classes)
 
     def staged_predict_proba(self, X):
         """Predict class probabilities for X.
@@ -819,30 +824,9 @@ def staged_predict_proba(self, X):
         X = self._validate_data(X)
 
         n_classes = self.n_classes_
-        proba = None
-        norm = 0.
-
-        for weight, estimator in zip(self.estimator_weights_,
-                                     self.estimators_):
-            norm += weight
-
-            if self.algorithm == 'SAMME.R':
-                # The weights are all 1. for SAMME.R
-                current_proba = _samme_proba(estimator, n_classes, X)
-            else:  # elif self.algorithm == "SAMME":
-                current_proba = estimator.predict_proba(X) * weight
-
-            if proba is None:
-                proba = current_proba
-            else:
-                proba += current_proba
-
-            real_proba = np.exp((1. / (n_classes - 1)) * (proba / norm))
-            normalizer = real_proba.sum(axis=1)[:, np.newaxis]
-            normalizer[normalizer == 0.0] = 1.0
-            real_proba /= normalizer
 
-            yield real_proba
+        for decision in self.staged_decision_function(X):
+            yield self._compute_proba_from_decision(decision, n_classes)
 
     def predict_log_proba(self, X):
         """Predict class log-probabilities for X.


Code review comments:
diff_hunk: @@ -748,6 +749,24 @@ class in ``classes_``, respectively.
             else:
                 yield pred / norm
 
+    @staticmethod
+    def _compute_proba_from_decision(decision, n_classes):
+        """Compute probabilities from the decision function.
+
+        This is based eq. (4) of [1] where:
+            p(y=c|X) = exp((1 / K-1) f_c(X)) / sum_k(exp((1 / K-1) f_k(X)))
+                     = softmax((1 / K-1) * f(X))
+
+        References
+        ----------
+        .. [1] J. Zhu, H. Zou, S. Rosset, T. Hastie, "Multi-class AdaBoost",
+               2009.
+        """
+        if n_classes == 2:
+            decision = np.vstack([-decision, decision]).T
+        decision *= 1 / (n_classes - 1)
path: sklearn/ensemble/weight_boosting.py
text: @user1:
use /= ?

diff_hunk: @@ -748,6 +749,24 @@ class in ``classes_``, respectively.
             else:
                 yield pred / norm
 
+    @staticmethod
+    def _compute_proba_from_decision(decision, n_classes):
+        """Compute probabilities from the decision function.
+
+        This is based eq. (4) of [1] where:
+            p(y=c|X) = exp((1 / K-1) f_c(X)) / sum_k(exp((1 / K-1) f_k(X)))
+                     = softmax((1 / K-1) * f(X))
+
+        References
+        ----------
+        .. [1] J. Zhu, H. Zou, S. Rosset, T. Hastie, "Multi-class AdaBoost",
+               2009.
+        """
+        if n_classes == 2:
+            decision = np.vstack([-decision, decision]).T
path: sklearn/ensemble/weight_boosting.py
text: @user1:
can you explain this? is it in the paper too?

@author:
My original thought was that we were keeping one of the 2 columns when computing the decision function. Actually, we do something a bit different:

```python
        if n_classes == 2:
            pred[:, 0] *= -1
            return pred.sum(axis=1)
```

You originally have symmetry between both classes. However, I should divide decision by 2 since that we are summing.