Allow MapieRegressor to use the optimal estimation strategy for the bounds of the prediction intervals

HISTORY.rst

@@ -5,6 +5,7 @@ History
 ##### (##########)
 ------------------
 * Integrate ConformityScore into MapieTimeSeriesRegressor.
+* Add (extend) the optimal estimation strategy for the bounds of the prediction intervals for regression via ConformityScore.
 * Add new checks for metrics calculations.
 * Fix reference for residual normalised score in documentation.
 

examples/regression/2-advanced-analysis/plot_timeseries_enbpi.py

@@ -36,6 +36,7 @@
 from sklearn.model_selection import RandomizedSearchCV, TimeSeriesSplit
 
 from mapie._typing import NDArray
+from mapie.conformity_scores import AbsoluteConformityScore
 from mapie.metrics import (regression_coverage_score,
                            regression_mean_width_score)
 from mapie.subsample import BlockBootstrap
@@ -111,6 +112,7 @@
     model,
     method="enbpi",
     cv=cv_mapietimeseries,
+    conformity_score=AbsoluteConformityScore(sym=False),
     agg_function="mean",
     n_jobs=-1,
 )

examples/regression/4-tutorials/plot_ts-tutorial.py

@@ -51,6 +51,7 @@ class that block bootstraps the training set.
 from sklearn.ensemble import RandomForestRegressor
 from sklearn.model_selection import RandomizedSearchCV, TimeSeriesSplit
 
+from mapie.conformity_scores import AbsoluteConformityScore
 from mapie.metrics import (regression_coverage_score,
                            regression_mean_width_score)
 from mapie.subsample import BlockBootstrap
@@ -187,8 +188,10 @@ class that block bootstraps the training set.
 cv_mapiets = BlockBootstrap(
     n_resamplings=10, n_blocks=10, overlapping=False, random_state=59
 )
+cs = AbsoluteConformityScore(sym=False)
 mapie_enbpi = MapieTimeSeriesRegressor(
-    model, method="enbpi", cv=cv_mapiets, agg_function="mean", n_jobs=-1
+    model, method="enbpi", cv=cv_mapiets, conformity_score=cs,
+    agg_function="mean", n_jobs=-1
 )
 
 ##############################################################################
@@ -216,8 +219,6 @@ class that block bootstraps the training set.
 y_pred_pfit = np.zeros(y_pred_npfit.shape)
 y_pis_pfit = np.zeros(y_pis_npfit.shape)
 conformity_scores_pfit = []
-lower_quantiles_pfit = []
-higher_quantiles_pfit = []
 y_pred_pfit[:gap], y_pis_pfit[:gap, :, :] = mapie_enbpi.predict(
     X_test.iloc[:gap, :], alpha=alpha, ensemble=True, optimize_beta=True
 )
@@ -236,8 +237,6 @@ class that block bootstraps the training set.
         optimize_beta=True
     )
     conformity_scores_pfit.append(mapie_enbpi.conformity_scores_)
-    lower_quantiles_pfit.append(mapie_enbpi.lower_quantiles_)
-    higher_quantiles_pfit.append(mapie_enbpi.higher_quantiles_)
 coverage_pfit = regression_coverage_score(
     y_test, y_pis_pfit[:, 0, 0], y_pis_pfit[:, 1, 0]
 )

mapie/conformity_scores/conformity_scores.py

@@ -254,14 +254,72 @@ def get_quantile(
         ])
         return quantile
 
+    @staticmethod
+    def _beta_optimize(
+        alpha_np: NDArray,
+        upper_bounds: NDArray,
+        lower_bounds: NDArray,
+    ) -> NDArray:
+        """
+        Minimize the width of the PIs, for a given difference of quantiles.
+
+        Parameters
+        ----------
+        alpha_np: NDArray
+            The quantiles to compute.
+
+        upper_bounds: NDArray
+            The array of upper values.
+
+        lower_bounds: NDArray
+            The array of lower values.
+
+        Returns
+        -------
+        NDArray
+            Array of betas minimizing the differences
+            ``(1-alpa+beta)-quantile - beta-quantile``.
+        """
+        beta_np = np.full(
+            shape=(len(lower_bounds), len(alpha_np)),
+            fill_value=np.nan,
+            dtype=float,
+        )
+
+        for ind_alpha, _alpha in enumerate(alpha_np):
+            betas = np.linspace(
+                _alpha / (len(lower_bounds) + 1),
+                _alpha,
+                num=len(lower_bounds),
+                endpoint=True,
+            )
+            one_alpha_beta = np_nanquantile(
+                upper_bounds.astype(float),
+                1 - _alpha + betas,
+                axis=1,
+                method="higher",
+            )
+            beta = np_nanquantile(
+                lower_bounds.astype(float),
+                betas,
+                axis=1,
+                method="lower",
+            )
+            beta_np[:, ind_alpha] = betas[
+                np.argmin(one_alpha_beta - beta, axis=0)
+            ]
+
+        return beta_np
+
     def get_bounds(
         self,
         X: ArrayLike,
         estimator: EnsembleEstimator,
         conformity_scores: NDArray,
         alpha_np: NDArray,
-        ensemble: bool,
-        method: str
+        ensemble: bool = False,
+        method: str = 'base',
+        optimize_beta: bool = False,
     ) -> Tuple[NDArray, NDArray, NDArray]:
         """
         Compute bounds of the prediction intervals from the observed values,
@@ -285,13 +343,22 @@ def get_bounds(
         ensemble: bool
             Boolean determining whether the predictions are ensembled or not.
 
+            By default ``False``.
+
         method: str
             Method to choose for prediction interval estimates.
             The ``"plus"`` method implies that the quantile is calculated
             after estimating the bounds, whereas the other methods
             (among the ``"naive"``, ``"base"`` or ``"minmax"`` methods,
             for example) do the opposite.
 
+            By default ``base``.
+
+        optimize_beta: bool
+            Whether to optimize the PIs' width or not.
+
+            By default ``False``.
+
         Returns
         -------
         Tuple[NDArray, NDArray, NDArray]
@@ -300,13 +367,33 @@ def get_bounds(
             (n_samples, n_alpha).
             - The upper bounds of the prediction intervals of shape
             (n_samples, n_alpha).
+
+        Raises
+        ------
+        ValueError
+            If beta optimisation with symmetrical conformity score function.
         """
+        if self.sym and optimize_beta:
+            raise ValueError(
+                "Beta optimisation cannot be used with "
+                + "symmetrical conformity score function."
+            )
+
         y_pred, y_pred_low, y_pred_up = estimator.predict(X, ensemble)
         signed = -1 if self.sym else 1
 
+        if optimize_beta:
+            beta_np = self._beta_optimize(
+                alpha_np,
+                conformity_scores.reshape(1, -1),
+                conformity_scores.reshape(1, -1),
+            )
+        else:
+            beta_np = alpha_np / 2
+
         if method == "plus":
-            alpha_low = alpha_np if self.sym else alpha_np / 2
-            alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np / 2
+            alpha_low = alpha_np if self.sym else beta_np
+            alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np + beta_np
 
             conformity_scores_low = self.get_estimation_distribution(
                 X, y_pred_low, signed * conformity_scores
@@ -322,8 +409,8 @@ def get_bounds(
             )
         else:
             quantile_search = "higher" if self.sym else "lower"
-            alpha_low = 1 - alpha_np if self.sym else alpha_np / 2
-            alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np / 2
+            alpha_low = 1 - alpha_np if self.sym else beta_np
+            alpha_up = 1 - alpha_np if self.sym else 1 - alpha_np + beta_np
 
             quantile_low = self.get_quantile(
                 conformity_scores, alpha_low, axis=0, method=quantile_search

mapie/estimator/estimator.py

@@ -496,9 +496,12 @@ def predict(
             if self.method == "minmax":
                 y_pred_multi_low = np.min(y_pred_multi, axis=1, keepdims=True)
                 y_pred_multi_up = np.max(y_pred_multi, axis=1, keepdims=True)
-            else:
+            elif self.method == "plus":
                 y_pred_multi_low = y_pred_multi
                 y_pred_multi_up = y_pred_multi
+            else:
+                y_pred_multi_low = y_pred[:, np.newaxis]
+                y_pred_multi_up = y_pred[:, np.newaxis]
 
             if ensemble:
                 y_pred = aggregate_all(self.agg_function, y_pred_multi)

mapie/regression/regression.py

@@ -522,6 +522,7 @@ def predict(
         X: ArrayLike,
         ensemble: bool = False,
         alpha: Optional[Union[float, Iterable[float]]] = None,
+        optimize_beta: bool = False,
     ) -> Union[NDArray, Tuple[NDArray, NDArray]]:
         """
         Predict target on new samples with confidence intervals.
@@ -561,6 +562,11 @@ def predict(
 
             By default ``None``.
 
+        optimize_beta: bool
+            Whether to optimize the PIs' width or not.
+
+            By default ``False``.
+
         Returns
         -------
         Union[NDArray, Tuple[NDArray, NDArray]]
@@ -582,6 +588,12 @@ def predict(
             return np.array(y_pred)
 
         else:
+            if optimize_beta and self.method != 'enbpi':
+                raise UserWarning(
+                    "Beta optimisation should only be used for "
+                    "the method='enbpi'."
+                )
+
             n = len(self.conformity_scores_)
             alpha_np = cast(NDArray, alpha)
             check_alpha_and_n_samples(alpha_np, n)
@@ -592,7 +604,8 @@ def predict(
                     self.estimator_,
                     self.conformity_scores_,
                     alpha_np,
-                    ensemble,
-                    self.method
+                    ensemble=ensemble,
+                    method=self.method,
+                    optimize_beta=optimize_beta
                 )
             return np.array(y_pred), np.stack([y_pred_low, y_pred_up], axis=1)