FEA add temperature scaling to CalibratedClassifierCV (scikit-learn#31068)

Co-authored-by: Christian Lorentzen <lorentzen.ch@gmail.com>

Author: virchan

doc/modules/calibration.rst

@@ -276,6 +276,35 @@ probabilities, the calibrated probabilities for each class
 are predicted separately. As those probabilities do not necessarily sum to
 one, a postprocessing is performed to normalize them.
 
+On the other hand, temperature scaling naturally supports multiclass
+predictions by working with logits and finally applying the softmax function.
+
+Temperature Scaling
+^^^^^^^^^^^^^^^^^^^
+
+For a multi-class classification problem with :math:`n` classes, temperature scaling
+[9]_, `method="temperature"`, produces class probabilities by modifying the softmax
+function with a temperature parameter :math:`T`:
+
+.. math::
+       \mathrm{softmax}\left(\frac{z}{T}\right) \,,
+
+where, for a given sample, :math:`z` is the vector of logits for each class as predicted
+by the estimator to be calibrated. In terms of scikit-learn's API, this corresponds to
+the output of :term:`decision_function` or to the logarithm of :term:`predict_proba`.
+Probabilities are converted to logits by first adding a tiny positive constant to avoid
+numerical issues with logarithm of zero, and then applying the natural logarithm.
+
+The parameter :math:`T` is learned by minimizing :func:`~sklearn.metrics.log_loss`,
+i.e. cross-entropy loss, on a hold-out (calibration) set. Note that :math:`T` does not
+affect the location of the maximum in the softmax output. Therefore, temperature scaling
+does not alter the accuracy of the calibrating estimator.
+
+The main advantage of temperature scaling over other calibration methods is that it
+provides a natural way to obtain (better) calibrated multi-class probabilities with
+just one free parameter in contrast to using a "One-vs-Rest" scheme that adds more
+parameters for each single class.
+
 .. rubric:: Examples
 
 * :ref:`sphx_glr_auto_examples_calibration_plot_calibration_curve.py`
@@ -324,3 +353,7 @@ one, a postprocessing is performed to normalize them.
        :doi:`"Statistical Foundations of Actuarial Learning and its Applications"
        <10.1007/978-3-031-12409-9>`
        Springer Actuarial
+
+.. [9] `On Calibration of Modern Neural Networks
+       <https://proceedings.mlr.press/v70/guo17a/guo17a.pdf>`_,
+       C. Guo, G. Pleiss, Y. Sun, & K. Q. Weinberger, ICML 2017.

doc/whats_new/upcoming_changes/sklearn.calibration/31068.feature.rst

@@ -0,0 +1,2 @@
+- Added temperature scaling in :class:`calibration.CalibratedClassifierCV`.
+  By :user:`Virgil Chan <virchan>`.

sklearn/_loss/loss.py

@@ -457,6 +457,20 @@ def constant_to_optimal_zero(self, y_true, sample_weight=None):
         """Calculate term dropped in loss.
 
         With this term added, the loss of perfect predictions is zero.
+
+        Parameters
+        ----------
+        y_true : array-like of shape (n_samples,)
+            Observed, true target values.
+
+        sample_weight : None or array of shape (n_samples,), default=None
+            Sample weights.
+
+        Returns
+        -------
+        constant : ndarray of shape (n_samples,)
+            Constant value to be added to raw predictions so that the loss
+            of perfect predictions becomes zero.
         """
         return np.zeros_like(y_true)
 
@@ -982,8 +996,16 @@ class HalfMultinomialLoss(BaseLoss):
     classes: If the full hessian for classes k and l and sample i is H_i_k_l,
     we calculate H_i_k_k, i.e. k=l.
 
-    Reference
-    ---------
+    Parameters
+    ----------
+    sample_weight : {None, ndarray}
+        If sample_weight is None, the hessian might be constant.
+
+    n_classes : {None, int}
+        The number of classes for classification, else None.
+
+    References
+    ----------
     .. [1] :arxiv:`Simon, Noah, J. Friedman and T. Hastie.
         "A Blockwise Descent Algorithm for Group-penalized Multiresponse and
         Multinomial Regression".
@@ -1015,6 +1037,19 @@ def fit_intercept_only(self, y_true, sample_weight=None):
 
         This is the softmax of the weighted average of the target, i.e. over
         the samples axis=0.
+
+        Parameters
+        ----------
+        y_true : array-like of shape (n_samples,)
+            Observed, true target values.
+
+        sample_weight : None or array of shape (n_samples,), default=None
+            Sample weights.
+
+        Returns
+        -------
+        raw_prediction : numpy scalar or array of shape (n_classes,)
+            Raw predictions of an intercept-only model.
         """
         out = np.zeros(self.n_classes, dtype=y_true.dtype)
         eps = np.finfo(y_true.dtype).eps