diff --git a/HISTORY.rst b/HISTORY.rst index 8c1de5085..9102a2351 100644 --- a/HISTORY.rst +++ b/HISTORY.rst @@ -5,6 +5,7 @@ History ##### (##########) ------------------ +* Add cross conformal option for quantile regression * Add split conformal option for regression and classification * Update check method for calibration diff --git a/mapie/quantile_regression.py b/mapie/quantile_regression.py index a87afd74f..6241c33e4 100644 --- a/mapie/quantile_regression.py +++ b/mapie/quantile_regression.py @@ -4,20 +4,24 @@ from typing import Iterable, List, Optional, Tuple, Union, cast import numpy as np +from joblib import Parallel, delayed from sklearn.base import RegressorMixin, clone from sklearn.linear_model import QuantileRegressor -from sklearn.model_selection import train_test_split +from sklearn.model_selection import BaseCrossValidator, ShuffleSplit from sklearn.pipeline import Pipeline -from sklearn.utils import check_random_state from sklearn.utils.validation import (_check_y, _num_samples, check_is_fitted, indexable) -from ._compatibility import np_quantile +from ._compatibility import np_nanquantile from ._typing import ArrayLike, NDArray +from .aggregation_functions import aggregate_all +from .conformity_scores import ConformityScore from .regression import MapieRegressor -from .utils import (check_alpha_and_n_samples, +from .utils import (check_alpha, check_alpha_and_n_samples, + check_conformity_score, check_cv, check_defined_variables_predict_cqr, check_estimator_fit_predict, check_lower_upper_bounds, + check_n_features_in, check_nan_in_aposteriori_prediction, check_null_weight, fit_estimator) @@ -25,8 +29,7 @@ class MapieQuantileRegressor(MapieRegressor): """ This class implements the conformalized quantile regression strategy as proposed by Romano et al. (2019) to make conformal predictions. - The only valid ``method`` is "quantile" and the only valid default - ``cv`` is "split". + The valid ``method`` and ``cv`` are the same as for MapieRegressor. Parameters ---------- @@ -35,14 +38,120 @@ class MapieQuantileRegressor(MapieRegressor): (i.e. with fit and predict methods), by default ``None``. If ``None``, estimator defaults to a ``QuantileRegressor`` instance. - method: str - Method to choose for prediction, in this case, the only valid method - is the "quantile" method. - - cv: Optional[str] - By default the value is set to None. In theory a split method is - implemented as it is needed to provided both a training and calibration - set. + method: str, optional + Method to choose for prediction interval estimates. + Choose among: + + - "naive", based on training set conformity scores, + - "base", based on validation sets conformity scores, + - "plus", based on validation conformity scores and + testing predictions, + - "minmax", based on validation conformity scores and + testing predictions (min/max among cross-validation clones). + + By default ``"plus"``. + + cv: Optional[Union[int, str, BaseCrossValidator]] + The cross-validation strategy for computing conformity scores. + It directly drives the distinction between jackknife and cv variants. + Choose among: + + - ``None``, to use the default 5-fold cross-validation + - integer, to specify the number of folds. + If equal to 1, does not involve cross-validation but a division + of the data into training and calibration subsets. The splitter + used is the following: ``sklearn.model_selection.ShuffleSplit``. + - CV splitter: any ``sklearn.model_selection.BaseCrossValidator`` + Main variants are: + - ``sklearn.model_selection.LeaveOneOut`` (jackknife), + - ``sklearn.model_selection.KFold`` (cross-validation), + - ``subsample.Subsample`` object (bootstrap). + - ``"split"``, does not involve cross-validation but a division + of the data into training and calibration subsets. The splitter + used is the following: ``sklearn.model_selection.ShuffleSplit``. + - ``"prefit"``, assumes that ``estimator`` has been fitted already, + and the ``method`` parameter is ignored. + All data provided in the ``fit`` method is then used + for computing conformity scores only. + At prediction time, quantiles of these conformity scores are used + to provide a prediction interval with fixed width. + The user has to take care manually that data for model fitting and + conformity scores estimate are disjoint. + + By default ``None``. + + test_size: Optional[Union[int, float]] + If float, should be between 0.0 and 1.0 and represent the proportion + of the dataset to include in the test split. If int, represents the + absolute number of test samples. If None, it will be set to 0.1. + + If cv is not ``"split"``, ``test_size`` is ignored. + + By default ``None``. + + n_jobs: Optional[int] + Number of jobs for parallel processing using joblib + via the "locky" backend. + If ``-1`` all CPUs are used. + If ``1`` is given, no parallel computing code is used at all, + which is useful for debugging. + For n_jobs below ``-1``, ``(n_cpus + 1 - n_jobs)`` are used. + None is a marker for `unset` that will be interpreted as ``n_jobs=1`` + (sequential execution). + + By default ``None``. + + agg_function : str + Determines how to aggregate predictions from perturbed models, both at + training and prediction time. + + If ``None``, it is ignored except if cv class is ``Subsample``, + in which case an error is raised. + If "mean" or "median", returns the mean or median of the predictions + computed from the out-of-folds models. + Note: if you plan to set the ``ensemble`` argument to ``True`` in the + ``predict`` method, you have to specify an aggregation function. + Otherwise an error would be raised. + + The Jackknife+ interval can be interpreted as an interval around the + median prediction, and is guaranteed to lie inside the interval, + unlike the single estimator predictions. + + When the cross-validation strategy is Subsample (i.e. for the + Jackknife+-after-Bootstrap method), this function is also used to + aggregate the training set in-sample predictions. + + If cv is ``"prefit"`` or ``"split"``, ``agg_function`` is ignored. + + By default ``"mean"``. + + verbose : int, optional + The verbosity level, used with joblib for multiprocessing. + The frequency of the messages increases with the verbosity level. + If it more than ``10``, all iterations are reported. + Above ``50``, the output is sent to stdout. + + By default ``0``. + + conformity_score : Optional[ConformityScore] + ConformityScore instance. + It defines the link between the observed values, the predicted ones + and the conformity scores. For instance, the default None value + correspondonds to a conformity score which assumes + y_obs = y_pred + conformity_score. + + - ``None``, to use the default ``AbsoluteConformityScore`` conformity + score + - ConformityScore: any ``ConformityScore`` class + + By default ``None``. + + random_state: Optional[Union[int, RandomState]] + Pseudo random number generator state used for random uniform sampling + for evaluation quantiles and prediction sets in cumulated_score. + Pass an int for reproducible output across multiple function calls. + + By default ``None``. alpha: float Between 0 and 1.0, represents the risk level of the confidence @@ -51,7 +160,7 @@ class MapieQuantileRegressor(MapieRegressor): intervals. ``alpha`` is the complement of the target coverage level. - By default 0.1. + By default ``0.1``. Attributes ---------- @@ -74,9 +183,6 @@ class MapieQuantileRegressor(MapieRegressor): quantile of 1 - alpha/2 - [:, 2]: maximum of those first two scores - n_calib_samples: int - Number of samples in the calibration dataset. - References ---------- Yaniv Romano, Evan Patterson and Emmanuel J. Candès. @@ -87,32 +193,35 @@ class MapieQuantileRegressor(MapieRegressor): -------- >>> import numpy as np >>> from mapie.quantile_regression import MapieQuantileRegressor - >>> X_train = np.array([[0], [1], [2], [3], [4], [5]]) - >>> y_train = np.array([5, 7.5, 9.5, 10.5, 12.5, 15]) - >>> X_calib = np.array([[0], [1], [2], [3], [4], [5], [6], [7], [8], [9]]) - >>> y_calib = np.array([5, 7, 9, 4, 8, 1, 5, 7.5, 9.5, 12]) - >>> mapie_reg = MapieQuantileRegressor().fit( - ... X_train, - ... y_train, - ... X_calib=X_calib, - ... y_calib=y_calib - ... ) - >>> y_pred, y_pis = mapie_reg.predict(X_train) + >>> X = np.array([[0], [1], [2], [3], [4], [5], [6], [7], [8], [9], [10]]) + >>> y = np.array([5, 7.5, 9.5, 10.5, 12.5, 15, 16, 17.5, 18.5, 20, 21]) + >>> mapie_reg = MapieQuantileRegressor(alpha=0.25, random_state=42) + >>> mapie_reg = mapie_reg.fit(X, y) + >>> y_pred, y_pis = mapie_reg.predict(X, alpha=0.25) >>> print(y_pis[:, :, 0]) - [[-8.16666667 19. ] - [-6.33333333 20.83333333] - [-4.5 22.66666667] - [-2.66666667 24.5 ] - [-0.83333333 26.33333333] - [ 1. 28.16666667]] + [[ 5. 6.4 ] + [ 6.5 8. ] + [ 8. 9.6 ] + [ 9.5 11.2 ] + [11. 12.9 ] + [12.66666667 14.6 ] + [14.33333333 16.3 ] + [16. 18. ] + [17.66666667 19.7 ] + [19.33333333 21.4 ] + [21. 23.1 ]] >>> print(y_pred) - [ 5. 7. 9. 11. 13. 15.] + [ 5.92857143 7.5 9.07142857 10.64285714 12.21428571 13.78571429 + 15.35714286 16.92857143 18.5 20.07142857 21.64285714] """ - valid_methods_ = ["quantile"] fit_attributes = [ + "single_estimator_", + "single_estimator_alpha_", "estimators_", + "k_", "conformity_scores_", - "n_calib_samples", + "conformity_score_function_", + "n_features_in_", ] quantile_estimator_params = { @@ -143,15 +252,27 @@ def __init__( List[Union[RegressorMixin, Pipeline]] ] ] = None, - method: str = "quantile", - cv: Optional[str] = None, + method: str = "plus", alpha: float = 0.1, + cv: Optional[Union[int, str, BaseCrossValidator]] = None, + test_size: Optional[Union[int, float]] = None, + n_jobs: Optional[int] = None, + agg_function: Optional[str] = "mean", + verbose: int = 0, + conformity_score: Optional[ConformityScore] = None, + random_state: Optional[Union[int, np.random.RandomState]] = None, ) -> None: super().__init__( estimator=estimator, method=method, + cv=cv, + test_size=test_size, + n_jobs=n_jobs, + agg_function=agg_function, + verbose=verbose, + conformity_score=conformity_score, + random_state=random_state ) - self.cv = cv self.alpha = alpha def _check_alpha( @@ -184,6 +305,7 @@ def _check_alpha( ------ ValueError If alpha is not a float. + ValueError If the value of alpha is not between 0 and 1.0. """ @@ -234,17 +356,21 @@ def _check_estimator( ------ ValueError If the estimator fit or predict methods. + ValueError We check if it's a known estimator that does quantile regression according to the dictionnary set quantile_estimator_params. This dictionnary will need to be updated with the latest new available estimators. + ValueError The estimator does not have the "loss_name" in its parameters and therefore can not be used as an estimator. + ValueError There is no quantile "loss_name" and therefore this estimator can not be used as a ``MapieQuantileRegressor``. + ValueError The parameter to set the alpha value does not exist in this estimator and therefore we cannot use it. @@ -254,6 +380,11 @@ def _check_estimator( solver="highs-ds", alpha=0.0, ) + + if check_cv(self.cv) == "prefit": + self._check_prefit_params(estimator) + return estimator + check_estimator_fit_predict(estimator) if isinstance(estimator, Pipeline): self._check_estimator(estimator[-1]) @@ -294,116 +425,11 @@ def _check_estimator( "The base model does not seem to be accepted" + " by MapieQuantileRegressor. \n" "Give a base model among: \n" - "``quantile_estimator_params.keys()``" + f"{self.quantile_estimator_params.keys()}" "Or, add your base model to" + " ``quantile_estimator_params``." ) - def _check_cv( - self, - cv: Optional[str] = None - ) -> str: - """ - Check if cv argument is None, "split" or "prefit". - - Parameters - ---------- - cv : Optional[str], optional - cv to check, by default ``None``. - - Returns - ------- - str - cv itself or a default "split". - - Raises - ------ - ValueError - Raises an error if the cv is anything else but the method "split" - or "prefit. - Only the split method has been implemented. - """ - if cv is None: - return "split" - if cv in ("split", "prefit"): - return cv - else: - raise ValueError( - "Invalid cv method, only valid method is ``split``." - ) - - def _check_calib_set( - self, - X: ArrayLike, - y: ArrayLike, - sample_weight: Optional[ArrayLike] = None, - X_calib: Optional[ArrayLike] = None, - y_calib: Optional[ArrayLike] = None, - calib_size: Optional[float] = 0.3, - random_state: Optional[Union[int, np.random.RandomState, None]] = None, - shuffle: Optional[bool] = True, - stratify: Optional[ArrayLike] = None, - ) -> Tuple[ - ArrayLike, ArrayLike, ArrayLike, ArrayLike, Optional[ArrayLike] - ]: - """ - Check if a calibration set has already been defined, if not, then - we define one using the `train_test_split` method. - - Parameters - ---------- - Same definition of parameters as for the ``fit`` method. - - Returns - ------- - Tuple[ArrayLike, ArrayLike, ArrayLike, ArrayLike, ArrayLike] - - [0]: ArrayLike of shape (n_samples_*(1-calib_size), n_features) - X_train - - [1]: ArrayLike of shape (n_samples_*(1-calib_size),) - y_train - - [2]: ArrayLike of shape (n_samples_*calib_size, n_features) - X_calib - - [3]: ArrayLike of shape (n_samples_*calib_size,) - y_calib - - [4]: ArrayLike of shape (n_samples_,) - sample_weight_train - - """ - if X_calib is None or y_calib is None: - if sample_weight is None: - X_train, X_calib, y_train, y_calib = train_test_split( - X, - y, - test_size=calib_size, - random_state=random_state, - shuffle=shuffle, - stratify=stratify - ) - sample_weight_train = sample_weight - else: - ( - X_train, - X_calib, - y_train, - y_calib, - sample_weight_train, - _, - ) = train_test_split( - X, - y, - sample_weight, - test_size=calib_size, - random_state=random_state, - shuffle=shuffle, - stratify=stratify - ) - else: - X_train, y_train, sample_weight_train = X, y, sample_weight - X_train, X_calib = cast(ArrayLike, X_train), cast(ArrayLike, X_calib) - y_train, y_calib = cast(ArrayLike, y_train), cast(ArrayLike, y_calib) - sample_weight_train = cast(ArrayLike, sample_weight_train) - return X_train, y_train, X_calib, y_calib, sample_weight_train - def _check_prefit_params( self, estimator: List[Union[RegressorMixin, Pipeline]], @@ -417,25 +443,21 @@ def _check_prefit_params( estimator : List[Union[RegressorMixin, Pipeline]] List of three prefitted estimators that should have pre-defined quantile levels of alpha/2, 1 - alpha/2 and 0.5. - X : ArrayLike of shape (n_samples, n_features) - Training data. - y : ArrayLike of shape (n_samples,) - Training labels. - X_calib : Optional[ArrayLike] of shape (n_calib_samples, n_features) - Calibration data. - y_calib : Optional[ArrayLike] of shape (n_calib_samples,) - Calibration labels. Raises ------ ValueError If a non-iterable variable is provided for estimator. + ValueError If less or more than three models are defined. + Warning If X and y are defined, then warning that they are not used. + ValueError If the calibration set is not defined. + Warning If the alpha is defined, warns the user that it must be set accordingly with the prefit estimators. @@ -455,17 +477,43 @@ def _check_prefit_params( " order [alpha/2, 1 - alpha/2, 0.5]." ) + def _pred_multi_alpha(self, X: ArrayLike) -> NDArray: + """ + TODO + Return a prediction per train sample for each test sample, by + aggregation with matrix ``k_``. + + Parameters + ---------- + X: NDArray of shape (n_samples_test, n_features) + Input data + + Returns + ------- + NDArray of shape (n_samples_test, n_samples_train) + """ + y_pred_alpha_list = [] + for i, est_list in enumerate(self.estimators_): + y_pred_multi = np.column_stack( + [e.predict(X) for e in est_list] + ) + # At this point, y_pred_multi is of shape + # (n_samples_test, n_estimators_). The method + # ``_aggregate_with_mask`` fits it to the right size + # thanks to the shape of k_. + y_pred_alpha_list.append( + self._aggregate_with_mask(y_pred_multi, self.k_[i]) + ) + + y_pred_alpha = np.array(y_pred_alpha_list) + + return y_pred_alpha + def fit( self, X: ArrayLike, y: ArrayLike, sample_weight: Optional[ArrayLike] = None, - X_calib: Optional[ArrayLike] = None, - y_calib: Optional[ArrayLike] = None, - calib_size: Optional[float] = 0.3, - random_state: Optional[Union[int, np.random.RandomState, None]] = None, - shuffle: Optional[bool] = True, - stratify: Optional[ArrayLike] = None, ) -> MapieQuantileRegressor: """ Fit estimator and compute residuals used for prediction intervals. @@ -479,8 +527,10 @@ def fit( ---------- X : ArrayLike of shape (n_samples, n_features) Training data. + y : ArrayLike of shape (n_samples,) Training labels. + sample_weight : Optional[ArrayLike] of shape (n_samples,) Sample weights for fitting the out-of-fold models. If None, then samples are equally weighted. @@ -490,124 +540,135 @@ def fit( If weights are non-uniform, residuals are still uniformly weighted. Note that the sample weight defined are only for the training, not for the calibration procedure. + By default ``None``. - X_calib : Optional[ArrayLike] of shape (n_calib_samples, n_features) - Calibration data. - y_calib : Optional[ArrayLike] of shape (n_calib_samples,) - Calibration labels. - calib_size : Optional[float] - If X_calib and y_calib are not defined, then the calibration - dataset is created with the split defined by calib_size. - random_state : int, RandomState instance or None, default=None - For the ``sklearn.model_selection.train_test_split`` documentation. - Controls the shuffling applied to the data before applying the - split. - Pass an int for reproducible output across multiple function calls. - See :term:`Glossary `. - shuffle : bool, default=True - For the ``sklearn.model_selection.train_test_split`` documentation. - Whether or not to shuffle the data before splitting. - If shuffle=False - then stratify must be None. - stratify : array-like, default=None - For the ``sklearn.model_selection.train_test_split`` documentation. - If not None, data is split in a stratified fashion, using this as - the class labels. - Read more in the :ref:`User Guide `. Returns ------- MapieQuantileRegressor The model itself. """ - self.cv = self._check_cv(cast(str, self.cv)) + # Checks + self._check_parameters() + cv = check_cv( + self.cv, test_size=self.test_size, random_state=self.random_state + ) + alpha = self._check_alpha(self.alpha) + estimator = self._check_estimator(self.estimator) + agg_function = self._check_agg_function(self.agg_function) + X, y = indexable(X, y) + y = _check_y(y) + sample_weight = cast(Optional[NDArray], sample_weight) + self.n_features_in_ = check_n_features_in(X, cv, estimator) + sample_weight, X, y = check_null_weight(sample_weight, X, y) + self.conformity_score_function_ = check_conformity_score( + self.conformity_score + ) + y = cast(NDArray, y) + n_samples = _num_samples(y) + check_alpha_and_n_samples(alpha, n_samples) + cv = cast(BaseCrossValidator, cv) # Initialization - self.estimators_: List[RegressorMixin] = [] - if self.cv == "prefit": - estimator = cast(List, self.estimator) - alpha = self._check_alpha(self.alpha) - self._check_prefit_params(estimator) - X_calib, y_calib = indexable(X, y) - - self.n_calib_samples = _num_samples(y_calib) - y_calib_preds = np.full( - shape=(3, self.n_calib_samples), - fill_value=np.nan - ) - for i, est in enumerate(estimator): - self.estimators_.append(est) - y_calib_preds[i] = est.predict(X_calib).ravel() - self.single_estimator_ = self.estimators_[2] - else: - # Checks - self._check_parameters() - checked_estimator = self._check_estimator(self.estimator) - alpha = self._check_alpha(self.alpha) - X, y = indexable(X, y) - random_state = check_random_state(random_state) - results = self._check_calib_set( - X, - y, - sample_weight, - X_calib, - y_calib, - calib_size, - random_state, - shuffle, - stratify, + self.single_estimator_alpha_: List[RegressorMixin] = [] + self.estimators_: List[List[RegressorMixin]] = [] + + def clone_estimator(_estimator, _alpha): + cloned_estimator_ = clone(_estimator) + if isinstance(_estimator, Pipeline): + alpha_name = self.quantile_estimator_params[ + _estimator[-1].__class__.__name__ + ]["alpha_name"] + _params = {alpha_name: _alpha} + cloned_estimator_[-1].set_params(**_params) + else: + alpha_name = self.quantile_estimator_params[ + _estimator.__class__.__name__ + ]["alpha_name"] + _params = {alpha_name: _alpha} + cloned_estimator_.set_params(**_params) + return cloned_estimator_ + + # Work + y_pred = np.full( + shape=(3, n_samples), + fill_value=np.nan, + dtype=float, + ) + self.conformity_scores_ = np.full( + shape=(3, n_samples), + fill_value=np.nan, + dtype=float, + ) + if self.cv != "prefit": + pred_matrix = np.full( + shape=(3, n_samples, cv.get_n_splits(X, y)), + fill_value=np.nan, + dtype=float, ) - X_train, y_train, X_calib, y_calib, sample_weight_train = results - X_train, y_train = indexable(X_train, y_train) - X_calib, y_calib = indexable(X_calib, y_calib) - y_train, y_calib = _check_y(y_train), _check_y(y_calib) - self.n_calib_samples = _num_samples(y_calib) - check_alpha_and_n_samples(self.alpha, self.n_calib_samples) - sample_weight_train, X_train, y_train = check_null_weight( - sample_weight_train, - X_train, - y_train + self.k_ = np.full( + shape=(3, n_samples, cv.get_n_splits(X, y)), + fill_value=np.nan, + dtype=float, ) - y_train = cast(NDArray, y_train) - - y_calib_preds = np.full( - shape=(3, self.n_calib_samples), - fill_value=np.nan + else: + self.k_ = np.full( + shape=(3, n_samples, 1), fill_value=np.nan, dtype=float ) - if isinstance(checked_estimator, Pipeline): - estimator = checked_estimator[-1] + for i, alpha_ in enumerate(alpha): + if self.cv == "prefit": + self.single_estimator_alpha_.append(estimator[i]) + y_pred[i] = self.single_estimator_alpha_[-1].predict(X) else: - estimator = checked_estimator - name_estimator = estimator.__class__.__name__ - alpha_name = self.quantile_estimator_params[ - name_estimator - ]["alpha_name"] - for i, alpha_ in enumerate(alpha): - cloned_estimator_ = clone(checked_estimator) - params = {alpha_name: alpha_} - if isinstance(checked_estimator, Pipeline): - cloned_estimator_[-1].set_params(**params) - else: - cloned_estimator_.set_params(**params) - self.estimators_.append(fit_estimator( - cloned_estimator_, X_train, y_train, sample_weight_train + cloned_estimator_ = clone_estimator(estimator, alpha_) + self.single_estimator_alpha_.append(fit_estimator( + cloned_estimator_, X, y, sample_weight )) - y_calib_preds[i] = self.estimators_[-1].predict(X_calib) - self.single_estimator_ = self.estimators_[2] - self.conformity_scores_ = np.full( - shape=(3, self.n_calib_samples), - fill_value=np.nan + if self.method == "naive": + y_pred[i] = self.single_estimator_alpha_[-1].predict(X) + else: + outputs = Parallel(n_jobs=self.n_jobs, + verbose=self.verbose)( + delayed(self._fit_and_predict_oof_model)( + clone_estimator(estimator, alpha_), + X, + y, + train_index, + val_index, + sample_weight, + ) + for train_index, val_index in cv.split(X) + ) + new_estimators, predictions, val_indices = map( + list, zip(*outputs) + ) + + self.estimators_.append(new_estimators) + + for j, val_ind in enumerate(val_indices): + pred_matrix[i, val_ind, j] = np.array( + predictions[j], dtype=float + ) + self.k_[i, val_ind, j] = 1 + check_nan_in_aposteriori_prediction(pred_matrix[i]) + + y_pred[i] = aggregate_all(agg_function, pred_matrix[i]) + + self.conformity_score_function_.sym = False + self.conformity_scores_[i] = ( + self.conformity_score_function_.get_conformity_scores( + y, + y_pred[i] + ) ) - self.conformity_scores_[0] = y_calib_preds[0] - y_calib - self.conformity_scores_[1] = y_calib - y_calib_preds[1] - self.conformity_scores_[2] = np.max( - [ - self.conformity_scores_[0], - self.conformity_scores_[1] - ], axis=0 - ) + + if isinstance(cv, ShuffleSplit): + self.single_estimator_ = self.estimators_[2][0] + else: + self.single_estimator_ = self.single_estimator_alpha_[2] + return self def predict( @@ -615,30 +676,44 @@ def predict( X: ArrayLike, ensemble: bool = False, alpha: Optional[Union[float, Iterable[float]]] = None, - symmetry: Optional[bool] = True, ) -> Union[NDArray, Tuple[NDArray, NDArray]]: """ Predict target on new samples with confidence intervals. - Residuals from the training set and predictions from the model clones - are central to the computation. - Prediction Intervals for a given ``alpha`` are deduced from the - quantile regression at the alpha values: alpha/2, 1 - (alpha/2) - while adding a constant based uppon their residuals. + Conformity scores from the training set and predictions + from the model clones are central to the computation. + Prediction Intervals for a given ``alpha`` are deduced from either + + - quantiles of conformity scores (naive and base methods), + - quantiles of (predictions +/- conformity scores) (plus method), + - quantiles of (max/min(predictions) +/- conformity scores) + (minmax method). Parameters ---------- - X : ArrayLike of shape (n_samples, n_features) + X: ArrayLike of shape (n_samples, n_features) Test data. - ensemble : bool - Ensemble has not been defined in predict and therefore should - will not have any effects in this method. - alpha : Optional[Union[float, Iterable[float]]] - For ``MapieQuantileRegresor`` the alpha has to be defined - directly in initial arguments of the class. - symmetry : Optional[bool], optional - Deciding factor to whether to find the quantile value for - each residuals separatly or to use the maximum of the two - combined. + + ensemble: bool + Boolean determining whether the predictions are ensembled or not. + If False, predictions are those of the model trained on the whole + training set. + If True, predictions from perturbed models are aggregated by + the aggregation function specified in the ``agg_function`` + attribute. + + If cv is ``"prefit"`` or ``"split"``, ``ensemble`` is ignored. + + By default ``False``. + + alpha: Optional[Union[float, Iterable[float]]] + Can be a float, a list of floats, or a ``ArrayLike`` of floats. + Between 0 and 1, represents the uncertainty of the confidence + interval. + Lower ``alpha`` produce larger (more conservative) prediction + intervals. + ``alpha`` is the complement of the target coverage level. + + By default ``None``. Returns ------- @@ -652,40 +727,93 @@ def predict( - [:, 0, :]: Lower bound of the prediction interval. - [:, 1, :]: Upper bound of the prediction interval. """ + # Checks check_is_fitted(self, self.fit_attributes) check_defined_variables_predict_cqr(ensemble, alpha) - alpha = self.alpha if symmetry else self.alpha/2 - check_alpha_and_n_samples(alpha, self.n_calib_samples) + self._check_ensemble(ensemble) + + y_pred = self.single_estimator_.predict(X) - n = self.n_calib_samples - q = (1 - (alpha)) * (1 + (1 / n)) + if alpha is None: + return np.array(y_pred) - y_preds = np.full( + n_samples = len(self.conformity_scores_[-1]) + alpha_np = cast(NDArray, check_alpha(alpha)) + check_alpha_and_n_samples(alpha_np, n_samples) + + y_pred = np.full( shape=(3, _num_samples(X)), fill_value=np.nan, dtype=float, ) - for i, est in enumerate(self.estimators_): - y_preds[i] = est.predict(X) - if symmetry: - quantile = np.full( - 2, - np_quantile( - self.conformity_scores_[2], q, method="higher" - ) - ) + for i, est in enumerate(self.single_estimator_alpha_): + y_pred[i] = est.predict(X) + + if self.method in self.no_agg_methods_ \ + or self.cv in self.no_agg_cv_: + y_pred_multi_low = y_pred[0, :, np.newaxis] + y_pred_multi_up = y_pred[1, :, np.newaxis] else: - quantile = np.array( - [ - np_quantile( - self.conformity_scores_[0], q, method="higher" - ), - np_quantile( - self.conformity_scores_[1], q, method="higher" + y_pred_multi = self._pred_multi_alpha(X) + + if self.method == "minmax": + y_pred_multi_low = np.min( + y_pred_multi[0], axis=1, keepdims=True + ) + y_pred_multi_up = np.max( + y_pred_multi[1], axis=1, keepdims=True + ) + else: + y_pred_multi_low = y_pred_multi[0] + y_pred_multi_up = y_pred_multi[1] + + if ensemble: + for i, est in enumerate(self.single_estimator_alpha_): + y_pred[i] = aggregate_all( + self.agg_function, y_pred_multi[i] ) - ] + + # compute distributions of lower and upper bounds + conformity_scores_low = self.conformity_scores_[0] + conformity_scores_up = self.conformity_scores_[1] + alpha_np = alpha_np/2 + + lower_bounds = ( + self.conformity_score_function_.get_estimation_distribution( + y_pred_multi_low, conformity_scores_low + ) + ) + upper_bounds = ( + self.conformity_score_function_.get_estimation_distribution( + y_pred_multi_up, conformity_scores_up ) - y_pred_low = y_preds[0][:, np.newaxis] - quantile[0] - y_pred_up = y_preds[1][:, np.newaxis] + quantile[1] - check_lower_upper_bounds(y_preds, y_pred_low, y_pred_up) - return y_preds[2], np.stack([y_pred_low, y_pred_up], axis=1) + ) + + # get desired confidence intervals according to alpha + y_pred_low = np.column_stack( + [ + np_nanquantile( + lower_bounds.astype(float), + _alpha, + axis=1, + method="lower", + ) + for _alpha in alpha_np + ] + ) + y_pred_up = np.column_stack( + [ + np_nanquantile( + upper_bounds.astype(float), + 1-_alpha, + axis=1, + method="higher", + ) + for _alpha in alpha_np + ] + ) + + for i in range(len(alpha_np)): + check_lower_upper_bounds(y_pred, y_pred_low[:, i], y_pred_up[:, i]) + + return y_pred[2], np.stack([y_pred_low, y_pred_up], axis=1) diff --git a/mapie/regression.py b/mapie/regression.py index a9b208f25..bcfc9910e 100644 --- a/mapie/regression.py +++ b/mapie/regression.py @@ -639,7 +639,7 @@ def predict( Parameters ---------- - X : ArrayLike of shape (n_samples, n_features) + X: ArrayLike of shape (n_samples, n_features) Test data. ensemble: bool @@ -680,16 +680,15 @@ def predict( # Checks check_is_fitted(self, self.fit_attributes) self._check_ensemble(ensemble) - alpha = cast(Optional[NDArray], check_alpha(alpha)) y_pred = self.single_estimator_.predict(X) - n = len(self.conformity_scores_) if alpha is None: return np.array(y_pred) - alpha_np = cast(NDArray, alpha) - check_alpha_and_n_samples(alpha_np, n) + n_samples = len(self.conformity_scores_) + alpha_np = cast(NDArray, check_alpha(alpha)) + check_alpha_and_n_samples(alpha_np, n_samples) if self.method in self.no_agg_methods_ \ or self.cv in self.no_agg_cv_: diff --git a/mapie/tests/test_quantile_regression.py b/mapie/tests/test_quantile_regression.py index d16226eda..e220291d5 100644 --- a/mapie/tests/test_quantile_regression.py +++ b/mapie/tests/test_quantile_regression.py @@ -1,24 +1,17 @@ from __future__ import annotations from inspect import signature -from typing import Any, Tuple +from typing import Any, Optional import numpy as np -import pandas as pd import pytest from sklearn.base import BaseEstimator, RegressorMixin, clone -from sklearn.compose import ColumnTransformer from sklearn.datasets import make_regression from sklearn.ensemble import GradientBoostingRegressor, RandomForestClassifier -from sklearn.impute import SimpleImputer from sklearn.linear_model import LinearRegression, QuantileRegressor -from sklearn.model_selection import KFold, LeaveOneOut, train_test_split -from sklearn.pipeline import Pipeline, make_pipeline -from sklearn.preprocessing import OneHotEncoder -from sklearn.utils.validation import check_is_fitted +from sklearn.model_selection import train_test_split from typing_extensions import TypedDict -from mapie._typing import NDArray from mapie.metrics import regression_coverage_score from mapie.quantile_regression import MapieQuantileRegressor @@ -33,25 +26,19 @@ random_state = 1 -X_train_toy, X_calib_toy, y_train_toy, y_calib_toy = train_test_split( - X_toy, - y_toy, - test_size=0.5, - random_state=random_state -) - qt = QuantileRegressor(solver="highs-ds") gb = GradientBoostingRegressor( - loss="quantile", - random_state=random_state - ) + loss="quantile", + random_state=random_state +) X, y = make_regression( n_samples=500, n_features=10, noise=1.0, random_state=random_state - ) +) + X_train, X_calib, y_train, y_calib = train_test_split( X, y, @@ -59,7 +46,6 @@ random_state=random_state ) -SYMMETRY = [True, False] ESTIMATOR = [qt, gb] Params = TypedDict( @@ -67,28 +53,34 @@ { "method": str, "alpha": float, + "random_state": Optional[int], }, ) STRATEGIES = { - "quantile_alpha2": Params(method="quantile", alpha=0.2), - "quantile_alpha3": Params(method="quantile", alpha=0.3), - "quantile_alpha4": Params(method="quantile", alpha=0.4), - "quantile_alpha8": Params(method="quantile", alpha=0.8), - } + "quantile_alpha2": + Params(method="plus", alpha=0.2, random_state=random_state), + "quantile_alpha3": + Params(method="plus", alpha=0.3, random_state=random_state), + "quantile_alpha4": + Params(method="plus", alpha=0.4, random_state=random_state), + "quantile_alpha8": + Params(method="plus", alpha=0.8, random_state=random_state), +} WIDTHS = { - "quantile_alpha2": 2.7360884795455576, - "quantile_alpha3": 2.185652142101473, - "quantile_alpha4": 1.731718678152845, - "quantile_alpha8": 0.66909752420949, + "quantile_alpha2": 2.744639410027496, + "quantile_alpha3": 2.100088918692794, + "quantile_alpha4": 1.7101104296999252, + "quantile_alpha8": 0.4907343583658844 } + COVERAGES = { - "quantile_alpha2": 0.834, - "quantile_alpha3": 0.738, - "quantile_alpha4": 0.646, - "quantile_alpha8": 0.264, + "quantile_alpha2": 0.84, + "quantile_alpha3": 0.732, + "quantile_alpha4": 0.608, + "quantile_alpha8": 0.206 } @@ -123,7 +115,7 @@ def predict(self, *args: Any) -> None: def test_default_parameters() -> None: """Test default values of input parameters.""" mapie_reg = MapieQuantileRegressor() - assert mapie_reg.method == "quantile" + assert mapie_reg.method == "plus" assert mapie_reg.cv is None assert mapie_reg.alpha == 0.1 @@ -140,13 +132,8 @@ def test_default_sample_weight() -> None: def test_default_parameters_estimator() -> None: """Test default values of estimator.""" mapie_reg = MapieQuantileRegressor() - mapie_reg.fit( - X_train, - y_train, - X_calib=X_calib, - y_calib=y_calib - ) - for estimator in mapie_reg.estimators_: + mapie_reg.fit(X, y) + for estimator in mapie_reg.single_estimator_alpha_: assert isinstance(estimator, QuantileRegressor) assert estimator.__dict__["solver"] == "highs-ds" @@ -159,13 +146,8 @@ def test_no_predict_fit_estimator() -> None: ): mapie_reg = MapieQuantileRegressor( estimator=NotFitPredictEstimator(alpha=0.2) - ) - mapie_reg.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) + ) + mapie_reg.fit(X_toy, y_toy) def test_no_para_loss_estimator() -> None: @@ -181,15 +163,8 @@ def test_no_para_loss_estimator() -> None: "loss_name": "noloss", "alpha_name": "alpha" } - mapie_reg.estimator = NoLossPamameterEstimator( - alpha=0.2 - ) - mapie_reg.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) + mapie_reg.estimator = NoLossPamameterEstimator(alpha=0.2) + mapie_reg.fit(X_toy, y_toy) def test_no_para_alpha_estimator() -> None: @@ -208,140 +183,26 @@ def test_no_para_alpha_estimator() -> None: mapie_reg.estimator = NoAlphaPamameterEstimator( alpha=0.2, loss="quantile" - ) - mapie_reg.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy ) + mapie_reg.fit(X_toy, y_toy) -@pytest.mark.parametrize("strategy", [*STRATEGIES]) @pytest.mark.parametrize("estimator", ESTIMATOR) -def test_valid_method(strategy: str, estimator: RegressorMixin) -> None: - """Test that valid strategies and estimators raise no error""" - mapie_reg = MapieQuantileRegressor( - estimator=estimator, - **STRATEGIES[strategy] - ) - mapie_reg.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) - check_is_fitted(mapie_reg, mapie_reg.fit_attributes) - assert mapie_reg.__dict__["method"] == "quantile" - - -@pytest.mark.parametrize("strategy", [*STRATEGIES]) -@pytest.mark.parametrize("estimator", ESTIMATOR) -@pytest.mark.parametrize("dataset", [ - (X_train, X_calib, y_train, y_calib), - (X_train_toy, X_calib_toy, y_train_toy, y_calib_toy) - ] -) -@pytest.mark.parametrize("symmetry", SYMMETRY) -def test_predict_output_shape( - strategy: str, - estimator: RegressorMixin, - dataset: Tuple[NDArray, NDArray, NDArray, NDArray], - symmetry: bool -) -> None: - """Test predict output shape.""" - mapie_reg = MapieQuantileRegressor( - estimator=estimator, - **STRATEGIES[strategy] - ) - (X_t, X_c, y_t, y_c) = dataset - mapie_reg.fit(X_t, y_t, X_calib=X_c, y_calib=y_c) - y_pred, y_pis = mapie_reg.predict(X_t, symmetry=symmetry) - assert y_pred.shape == (X_t.shape[0],) - assert y_pis[:, 0, 0].shape == (X_t.shape[0],) - assert y_pis[:, 1, 0].shape == (X_t.shape[0],) - - -@pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_results_with_constant_sample_weights( - strategy: str, -) -> None: - """ - Test predictions when sample weights are None - or constant with different values. - """ - n_samples = len(X_train) - mapie0 = MapieQuantileRegressor( - estimator=qt, - **STRATEGIES[strategy] - ) - mapie1 = MapieQuantileRegressor( - estimator=qt, - **STRATEGIES[strategy] - ) - mapie2 = MapieQuantileRegressor( - estimator=qt, - **STRATEGIES[strategy] - ) - mapie0.fit( - X_train, - y_train, - X_calib=X_calib, - y_calib=y_calib, - sample_weight=None - ) - mapie1.fit( - X_train, - y_train, - X_calib=X_calib, - y_calib=y_calib, - sample_weight=np.ones(shape=n_samples) - ) - mapie2.fit( - X_train, - y_train, - X_calib=X_calib, - y_calib=y_calib, - sample_weight=np.ones(shape=n_samples) * 5 - ) - - np.testing.assert_allclose( - mapie0.conformity_scores_, - mapie1.conformity_scores_ - ) - np.testing.assert_allclose( - mapie0.conformity_scores_, - mapie2.conformity_scores_ - ) - - y_pred0, y_pis0 = mapie0.predict(X) - y_pred1, y_pis1 = mapie1.predict(X) - y_pred2, y_pis2 = mapie2.predict(X) - np.testing.assert_allclose(y_pred0, y_pred1) - np.testing.assert_allclose(y_pred1, y_pred2) - np.testing.assert_allclose(y_pis0, y_pis1) - np.testing.assert_allclose(y_pis1, y_pis2) - - -@pytest.mark.parametrize("estimator", ESTIMATOR) -@pytest.mark.parametrize("symmetry", SYMMETRY) -def test_results_for_same_alpha( - estimator: RegressorMixin, - symmetry: bool -) -> None: +def test_results_for_same_alpha(estimator: RegressorMixin) -> None: """ Test that predictions and intervals are similar with two equal values of alpha. """ mapie_reg = MapieQuantileRegressor( estimator=estimator, - alpha=0.2 + alpha=0.2, + random_state=random_state ) mapie_reg_clone = clone(mapie_reg) - mapie_reg.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - mapie_reg_clone.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - y_pred, y_pis = mapie_reg.predict(X, symmetry=symmetry) - y_pred_clone, y_pis_clone = mapie_reg_clone.predict(X, symmetry=symmetry) + mapie_reg.fit(X, y) + mapie_reg_clone.fit(X, y) + y_pred, y_pis = mapie_reg.predict(X, alpha=0.2) + y_pred_clone, y_pis_clone = mapie_reg_clone.predict(X, alpha=0.2) np.testing.assert_allclose(y_pred, y_pred_clone) np.testing.assert_allclose(y_pis[:, 0, 0], y_pis_clone[:, 0, 0]) np.testing.assert_allclose(y_pis[:, 1, 0], y_pis_clone[:, 1, 0]) @@ -355,7 +216,7 @@ def test_wrong_alphas_types(alphas: float) -> None: match=r".*Invalid alpha. Allowed values are float.*", ): mapie_reg = MapieQuantileRegressor(alpha=alphas) - mapie_reg.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) + mapie_reg.fit(X, y) @pytest.mark.parametrize("alphas", [1.0, 1.6, 1.95, 5.0, -0.1, -0.001, -10.0]) @@ -366,7 +227,7 @@ def test_wrong_alphas(alphas: float) -> None: match=r".*Invalid alpha. Allowed values are between 0 and .*", ): mapie_reg = MapieQuantileRegressor(alpha=alphas) - mapie_reg.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) + mapie_reg.fit(X, y) def test_estimators_quantile_function() -> None: @@ -378,67 +239,7 @@ def test_estimators_quantile_function() -> None: mapie_reg = MapieQuantileRegressor( estimator=GradientBoostingRegressor() ) - mapie_reg.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - - -@pytest.mark.parametrize("cv", [-1, 2, KFold(), LeaveOneOut()]) -def test_invalid_cv(cv: Any) -> None: - """Test that valid cv raise errors.""" - with pytest.raises( - ValueError, - match=r".*Invalid cv method.*", - ): - mapie = MapieQuantileRegressor(cv=cv) - mapie.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) - - -@pytest.mark.parametrize("cv", [None, "split"]) -def test_valid_cv(cv: Any) -> None: - """Test that valid cv raise no errors.""" - mapie = MapieQuantileRegressor(cv=cv) - mapie.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) - - -def test_calib_dataset_is_none() -> None: - """Test that the fit method works when X_calib or y_calib is None.""" - mapie = MapieQuantileRegressor() - mapie.fit(X, y, calib_size=0.5) - mapie.predict(X) - - -def test_calib_dataset_is_none_with_sample_weight() -> None: - """ - Test that the fit method works with calib dataset defined is None - with sample weights. - """ - mapie = MapieQuantileRegressor() - mapie.fit(X, y, calib_size=0.5, sample_weight=np.ones(X.shape[0])) - mapie.predict(X) - - -def test_calib_dataset_is_none_vs_defined() -> None: - """ - Test that for the same results whether you split before - or in the fit method. - """ - mapie = MapieQuantileRegressor() - mapie_defined = clone(mapie) - mapie.fit(X, y, calib_size=0.5, random_state=random_state) - mapie_defined.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - y_pred, y_pis = mapie.predict(X) - y_pred_defined, y_pis_defined = mapie_defined.predict(X) - np.testing.assert_allclose(y_pred, y_pred_defined, rtol=1e-2) - np.testing.assert_allclose(y_pis, y_pis_defined, rtol=1e-2) + mapie_reg.fit(X, y) @pytest.mark.parametrize("est", [RandomForestClassifier(), LinearRegression()]) @@ -452,12 +253,7 @@ def test_estimators_not_in_list(est: RegressorMixin) -> None: match=r".*The base model does not seem to be accepted by.*", ): mapie_reg = MapieQuantileRegressor(estimator=est) - mapie_reg.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) + mapie_reg.fit(X_toy, y_toy) def test_for_small_dataset() -> None: @@ -472,32 +268,8 @@ def test_for_small_dataset() -> None: ) mapie_reg.fit( np.array([1, 2, 3]), - np.array([2, 2, 3]), - X_calib=np.array([3, 5]), - y_calib=np.array([2, 3]) - ) - - -@pytest.mark.parametrize("strategy", [*STRATEGIES]) -@pytest.mark.parametrize("estimator", ESTIMATOR) -@pytest.mark.parametrize("dataset", [ - (X_train, X_calib, y_train, y_calib), - (X_train_toy, X_calib_toy, y_train_toy, y_calib_toy) -]) -def test_conformity_len( - strategy: str, - estimator: RegressorMixin, - dataset: Tuple[NDArray, NDArray, NDArray, NDArray], -) -> None: - """Test conformity scores output shape.""" - (X_t, X_c, y_t, y_c) = dataset - n_samples = int(len(X_c)) - mapie_regressor = MapieQuantileRegressor( - estimator=estimator, - **STRATEGIES[strategy] + np.array([2, 2, 3]) ) - mapie_regressor.fit(X_t, y_t, X_calib=X_c, y_calib=y_c) - assert mapie_regressor.conformity_scores_[0].shape[0] == n_samples # Working but want to add both symmetry and different estimators @@ -507,9 +279,10 @@ def test_linear_regression_results(strategy: str) -> None: Test expected prediction intervals for a different strategies. """ + alpha = STRATEGIES[strategy]['alpha'] mapie = MapieQuantileRegressor(**STRATEGIES[strategy]) - mapie.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - _, y_pis = mapie.predict(X) + mapie.fit(X, y) + _, y_pis = mapie.predict(X, alpha=alpha) y_pred_low, y_pred_up = y_pis[:, 0, 0], y_pis[:, 1, 0] width_mean = (y_pred_up - y_pred_low).mean() coverage = regression_coverage_score(y, y_pred_low, y_pred_up) @@ -517,6 +290,20 @@ def test_linear_regression_results(strategy: str) -> None: np.testing.assert_allclose(coverage, COVERAGES[strategy], rtol=1e-2) +def test_error_model_prefit() -> None: + """ + Check that the model is an iterable object for "prefit". + """ + with pytest.raises( + ValueError, + match=r".*Estimator for prefit must be an iterable object*" + ): + mapie_reg = MapieQuantileRegressor( + estimator=object(), cv="prefit" + ) + mapie_reg.fit(X, y) + + def test_quantile_prefit_three_estimators() -> None: """ Test that there is a list with three estimators provided for @@ -534,10 +321,7 @@ def test_quantile_prefit_three_estimators() -> None: estimator=list_estimators, cv="prefit" ) - mapie_reg.fit( - X_calib, - y_calib - ) + mapie_reg.fit(X_calib, y_calib) def test_prefit_no_fit_predict() -> None: @@ -558,10 +342,7 @@ def test_prefit_no_fit_predict() -> None: cv="prefit", alpha=0.3 ) - mapie_reg.fit( - X_calib, - y_calib - ) + mapie_reg.fit(X_calib, y_calib) def test_non_trained_estimator() -> None: @@ -581,10 +362,7 @@ def test_non_trained_estimator() -> None: cv="prefit", alpha=0.3 ) - mapie_reg.fit( - X_calib, - y_calib - ) + mapie_reg.fit(X_calib, y_calib) def test_warning_alpha_prefit() -> None: @@ -605,10 +383,7 @@ def test_warning_alpha_prefit() -> None: cv="prefit", alpha=0.3 ) - mapie_reg.fit( - X_calib, - y_calib - ) + mapie_reg.fit(X_calib, y_calib) @pytest.mark.parametrize("alpha", [0.05, 0.1, 0.2, 0.3]) @@ -623,22 +398,19 @@ def test_prefit_and_non_prefit_equal(alpha: float) -> None: est = clone(qt) params = {"quantile": alpha_} est.set_params(**params) - est.fit(X_train, y_train) + est.fit(X, y) list_estimators.append(est) mapie_reg_prefit = MapieQuantileRegressor( estimator=list_estimators, cv="prefit", alpha=alpha ) - mapie_reg_prefit.fit(X_calib, y_calib) - y_pred_prefit, y_pis_prefit = mapie_reg_prefit.predict(X) + mapie_reg_prefit.fit(X, y) + y_pred_prefit, y_pis_prefit = mapie_reg_prefit.predict(X, alpha=alpha) - mapie_reg = MapieQuantileRegressor( - estimator=qt, - alpha=alpha - ) - mapie_reg.fit(X_train, y_train, X_calib=X_calib, y_calib=y_calib) - y_pred, y_pis = mapie_reg.predict(X) + mapie_reg = MapieQuantileRegressor(estimator=qt, alpha=alpha) + mapie_reg.fit(X, y) + y_pred, y_pis = mapie_reg.predict(X, alpha=alpha) np.testing.assert_allclose(y_pred_prefit, y_pred) np.testing.assert_allclose(y_pis_prefit, y_pis) @@ -665,7 +437,8 @@ def test_prefit_different_type_list_tuple_array(alpha: float) -> None: alpha=alpha ) mapie_reg_prefit_list.fit(X_calib, y_calib) - y_pred_prefit_list, y_pis_prefit_list = mapie_reg_prefit_list.predict(X) + y_pred_prefit_list, y_pis_prefit_list = \ + mapie_reg_prefit_list.predict(X, alpha=alpha) mapie_reg_prefit_tuple = MapieQuantileRegressor( estimator=tuple(list_estimators), @@ -673,7 +446,8 @@ def test_prefit_different_type_list_tuple_array(alpha: float) -> None: alpha=alpha ) mapie_reg_prefit_tuple.fit(X_calib, y_calib) - y_pred_prefit_tuple, y_pis_prefit_tuple = mapie_reg_prefit_tuple.predict(X) + y_pred_prefit_tuple, y_pis_prefit_tuple = \ + mapie_reg_prefit_tuple.predict(X, alpha=alpha) mapie_reg_prefit_array = MapieQuantileRegressor( estimator=np.array(list_estimators), @@ -681,54 +455,11 @@ def test_prefit_different_type_list_tuple_array(alpha: float) -> None: alpha=alpha ) mapie_reg_prefit_array.fit(X_calib, y_calib) - y_pred_prefit_array, y_pis_prefit_array = mapie_reg_prefit_array.predict(X) + y_pred_prefit_array, y_pis_prefit_array = \ + mapie_reg_prefit_array.predict(X, alpha=alpha) np.testing.assert_allclose(y_pred_prefit_list, y_pred_prefit_tuple) np.testing.assert_allclose(y_pis_prefit_list, y_pis_prefit_tuple) np.testing.assert_allclose(y_pred_prefit_list, y_pred_prefit_array) np.testing.assert_allclose(y_pis_prefit_list, y_pis_prefit_array) - - -@pytest.mark.parametrize("estimator", ESTIMATOR) -def test_pipeline_compatibility(estimator: RegressorMixin) -> None: - """Check that MAPIE works on pipeline based on pandas dataframes""" - X = pd.DataFrame( - { - "x_cat": ["A", "A", "B", "A", "A", "B", "A", "B", "B", "B"], - "x_num": [0, 1, 1, 4, np.nan, 5, 4, 3, np.nan, 3], - "y": [5, 7, 3, 9, 10, 8, 9, 7, 9, 8] - } - ) - y = pd.Series([5, 7, 3, 9, 10, 8, 9, 7, 10, 5]) - X_train_toy, X_calib_toy, y_train_toy, y_calib_toy = train_test_split( - X, - y, - test_size=0.5, - random_state=random_state - ) - numeric_preprocessor = Pipeline( - [ - ("imputer", SimpleImputer(strategy="mean")), - ] - ) - categorical_preprocessor = Pipeline( - steps=[ - ("encoding", OneHotEncoder(handle_unknown="ignore")) - ] - ) - preprocessor = ColumnTransformer( - [ - ("cat", categorical_preprocessor, ["x_cat"]), - ("num", numeric_preprocessor, ["x_num"]) - ] - ) - pipe = make_pipeline(preprocessor, estimator) - mapie = MapieQuantileRegressor(pipe, alpha=0.4) - mapie.fit( - X_train_toy, - y_train_toy, - X_calib=X_calib_toy, - y_calib=y_calib_toy - ) - mapie.predict(X) diff --git a/mapie/tests/test_regression.py b/mapie/tests/test_regression.py index cda243c3e..886d1901c 100644 --- a/mapie/tests/test_regression.py +++ b/mapie/tests/test_regression.py @@ -8,30 +8,57 @@ import pytest from sklearn.compose import ColumnTransformer from sklearn.datasets import make_regression -from sklearn.dummy import DummyRegressor from sklearn.impute import SimpleImputer -from sklearn.linear_model import LinearRegression -from sklearn.model_selection import (KFold, LeaveOneOut, ShuffleSplit, - train_test_split) +from sklearn.base import BaseEstimator +from sklearn.linear_model import LinearRegression, QuantileRegressor +from sklearn.model_selection import KFold, LeaveOneOut, ShuffleSplit from sklearn.pipeline import Pipeline, make_pipeline from sklearn.preprocessing import OneHotEncoder from sklearn.utils.validation import check_is_fitted from typing_extensions import TypedDict -from mapie._typing import ArrayLike, NDArray -from mapie.aggregation_functions import aggregate_all +from mapie._typing import NDArray from mapie.conformity_scores import (AbsoluteConformityScore, ConformityScore, GammaConformityScore) from mapie.metrics import regression_coverage_score from mapie.regression import MapieRegressor +from mapie.quantile_regression import MapieQuantileRegressor from mapie.subsample import Subsample X_toy = np.array([0, 1, 2, 3, 4, 5]).reshape(-1, 1) y_toy = np.array([5, 7, 9, 11, 13, 15]) + +X_toy, y_toy = make_regression( + n_samples=50, n_features=10, noise=1.0, random_state=1 +) X, y = make_regression( n_samples=500, n_features=10, noise=1.0, random_state=1 ) -k = np.ones(shape=(5, X.shape[1])) + + +MAPIE_ESTIMATORS = [ + MapieRegressor, + MapieQuantileRegressor +] + +MAPIE_SINGLE_ESTIMATORS = [ + MapieRegressor +] + +BASE_ESTIMATORS = { + "MapieRegressor": LinearRegression, + "MapieQuantileRegressor": QuantileRegressor +} + +PREFIT_ESTIMATORS = { + "MapieRegressor": LinearRegression().fit(X, y), + "MapieQuantileRegressor": [ + QuantileRegressor(quantile=0.05).fit(X, y), + QuantileRegressor(quantile=0.50).fit(X, y), + QuantileRegressor(quantile=0.95).fit(X, y) + ] +} + METHODS = ["naive", "base", "plus", "minmax"] random_state = 1 @@ -46,6 +73,7 @@ "random_state": Optional[int], }, ) + STRATEGIES = { "naive": Params( method="naive", @@ -127,102 +155,149 @@ } WIDTHS = { - "naive": 3.81, - "split": 3.87, - "jackknife": 3.89, - "jackknife_plus": 3.90, - "jackknife_minmax": 3.96, - "cv": 3.85, - "cv_plus": 3.90, - "cv_minmax": 4.04, - "prefit": 4.81, - "cv_plus_median": 3.90, - "jackknife_plus_ab": 3.90, - "jackknife_minmax_ab": 4.13, - "jackknife_plus_median_ab": 3.87, -} - + "MapieRegressor": { + "naive": 3.81, + "split": 3.87, + "jackknife": 3.89, + "jackknife_plus": 3.9, + "jackknife_minmax": 3.96, + "cv": 3.84, + "cv_plus": 3.90, + "cv_minmax": 4.04, + "prefit": 4.81, + "cv_plus_median": 3.9, + "jackknife_plus_ab": 3.89, + "jackknife_minmax_ab": 4.14, + "jackknife_plus_median_ab": 3.88 + }, + "MapieQuantileRegressor": { + "naive": 3.56, + "split": 3.94, + "jackknife": 4.03, + "jackknife_plus": 3.96, + "jackknife_minmax": 4.27, + "cv": 4.29, + "cv_plus": 4.22, + "cv_minmax": 4.71, + "prefit": 4.81, + "cv_plus_median": 3.9, + "jackknife_plus_ab": 4.01, + "jackknife_minmax_ab": 4.56, + "jackknife_plus_median_ab": 4.1}} COVERAGES = { - "naive": 0.952, - "split": 0.952, - "jackknife": 0.952, - "jackknife_plus": 0.952, - "jackknife_minmax": 0.952, - "cv": 0.958, - "cv_plus": 0.956, - "cv_minmax": 0.966, - "prefit": 0.980, - "cv_plus_median": 0.954, - "jackknife_plus_ab": 0.952, - "jackknife_minmax_ab": 0.970, - "jackknife_plus_median_ab": 0.960, + "MapieRegressor": { + "naive": 0.952, + "split": 0.952, + "jackknife": 0.952, + "jackknife_plus": 0.952, + "jackknife_minmax": 0.952, + "cv": 0.952, + "cv_plus": 0.954, + "cv_minmax": 0.962, + "prefit": 0.98, + "cv_plus_median": 0.954, + "jackknife_plus_ab": 0.954, + "jackknife_minmax_ab": 0.968, + "jackknife_plus_median_ab": 0.952}, + "MapieQuantileRegressor": { + "naive": 0.952, + "split": 0.968, + "jackknife": 0.97, + "jackknife_plus": 0.966, + "jackknife_minmax": 0.974, + "cv": 0.972, + "cv_plus": 0.972, + "cv_minmax": 0.984, + "prefit": 0.98, + "cv_plus_median": 0.954, + "jackknife_plus_ab": 0.962, + "jackknife_minmax_ab": 0.984, + "jackknife_plus_median_ab": 0.972 + } } -def test_default_parameters() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_default_parameters(estimator: BaseEstimator) -> None: """Test default values of input parameters.""" - mapie_reg = MapieRegressor() - assert mapie_reg.agg_function == "mean" - assert mapie_reg.method == "plus" + mapie_est = estimator() + assert mapie_est.agg_function == "mean" + assert mapie_est.method == "plus" +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_valid_estimator(strategy: str) -> None: +def test_valid_estimator(estimator: BaseEstimator, strategy: str) -> None: """Test that valid estimators are not corrupted, for all strategies.""" - mapie_reg = MapieRegressor( - estimator=DummyRegressor(), **STRATEGIES[strategy] - ) + base_estimator = BASE_ESTIMATORS[estimator.__name__] + mapie_reg = estimator(estimator=base_estimator(), **STRATEGIES[strategy]) mapie_reg.fit(X_toy, y_toy) - assert isinstance(mapie_reg.single_estimator_, DummyRegressor) - for estimator in mapie_reg.estimators_: - assert isinstance(estimator, DummyRegressor) + assert isinstance(mapie_reg.single_estimator_, base_estimator) + for elt in mapie_reg.estimators_: + assert isinstance(elt, base_estimator) or \ + list(map(type, elt)).count(base_estimator) == len(elt) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("method", METHODS) -def test_valid_method(method: str) -> None: +def test_valid_method(estimator: BaseEstimator, method: str) -> None: """Test that valid methods raise no errors.""" - mapie_reg = MapieRegressor(method=method) + mapie_reg = estimator(method=method) mapie_reg.fit(X_toy, y_toy) check_is_fitted(mapie_reg, mapie_reg.fit_attributes) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("agg_function", ["dummy", 0, 1, 2.5, [1, 2]]) -def test_invalid_agg_function(agg_function: Any) -> None: +def test_invalid_agg_function( + estimator: BaseEstimator, agg_function: Any +) -> None: """Test that invalid agg_functions raise errors.""" - mapie_reg = MapieRegressor(agg_function=agg_function) + mapie_reg = estimator(agg_function=agg_function) with pytest.raises(ValueError, match=r".*Invalid aggregation function.*"): mapie_reg.fit(X_toy, y_toy) - mapie_reg = MapieRegressor(agg_function=None) + +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_invalid_agg_function_as_none(estimator: BaseEstimator) -> None: + """Test that invalid agg_functions raise errors.""" + mapie_reg = estimator(agg_function=None) with pytest.raises(ValueError, match=r".*If ensemble is True*"): mapie_reg.fit(X_toy, y_toy) mapie_reg.predict(X_toy, ensemble=True) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("agg_function", [None, "mean", "median"]) -def test_valid_agg_function(agg_function: str) -> None: +def test_valid_agg_function( + estimator: BaseEstimator, agg_function: str +) -> None: """Test that valid agg_functions raise no errors.""" - mapie_reg = MapieRegressor(agg_function=agg_function) + mapie_reg = estimator(agg_function=agg_function) mapie_reg.fit(X_toy, y_toy) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize( "cv", [None, -1, 2, KFold(), LeaveOneOut(), - ShuffleSplit(n_splits=1), "prefit", "split"] + ShuffleSplit(n_splits=1), "split", "prefit"] ) -def test_valid_cv(cv: Any) -> None: +def test_valid_cv(estimator: BaseEstimator, cv: Any) -> None: """Test that valid cv raise no errors.""" - model = LinearRegression() - model.fit(X_toy, y_toy) - mapie_reg = MapieRegressor(estimator=model, cv=cv) + if cv != "prefit": + model = BASE_ESTIMATORS[estimator.__name__]() + else: + model = PREFIT_ESTIMATORS[estimator.__name__] + mapie_reg = estimator(estimator=model, cv=cv) mapie_reg.fit(X_toy, y_toy) mapie_reg.predict(X_toy, alpha=0.5) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("cv", [100, 200, 300]) -def test_too_large_cv(cv: Any) -> None: +def test_too_large_cv(estimator: BaseEstimator, cv: Any) -> None: """Test that too large cv raise sklearn errors.""" - mapie_reg = MapieRegressor(cv=cv) + mapie_reg = estimator(cv=cv) with pytest.raises( ValueError, match=rf".*Cannot have number of splits n_splits={cv} greater.*", @@ -230,14 +305,16 @@ def test_too_large_cv(cv: Any) -> None: mapie_reg.fit(X_toy, y_toy) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -@pytest.mark.parametrize("dataset", [(X, y), (X_toy, y_toy)]) +@pytest.mark.parametrize("dataset", [(X_toy, y_toy)]) @pytest.mark.parametrize("alpha", [0.2, [0.2, 0.4], (0.2, 0.4)]) def test_predict_output_shape( - strategy: str, alpha: Any, dataset: Tuple[NDArray, NDArray] + estimator: BaseEstimator, strategy: str, alpha: Any, + dataset: Tuple[NDArray, NDArray] ) -> None: """Test predict output shape.""" - mapie_reg = MapieRegressor(**STRATEGIES[strategy]) + mapie_reg = estimator(**STRATEGIES[strategy]) (X, y) = dataset mapie_reg.fit(X, y) y_pred, y_pis = mapie_reg.predict(X, alpha=alpha) @@ -273,28 +350,32 @@ def test_same_results_prefit_split() -> None: np.testing.assert_allclose(y_pis_1[:, 1, 0], y_pis_2[:, 1, 0]) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_results_for_same_alpha(strategy: str) -> None: +def test_results_for_same_alpha( + estimator: BaseEstimator, strategy: str +) -> None: """ Test that predictions and intervals are similar with two equal values of alpha. """ - mapie_reg = MapieRegressor(**STRATEGIES[strategy]) + mapie_reg = estimator(**STRATEGIES[strategy]) mapie_reg.fit(X, y) _, y_pis = mapie_reg.predict(X, alpha=[0.1, 0.1]) np.testing.assert_allclose(y_pis[:, 0, 0], y_pis[:, 0, 1]) np.testing.assert_allclose(y_pis[:, 1, 0], y_pis[:, 1, 1]) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) @pytest.mark.parametrize( "alpha", [np.array([0.05, 0.1]), [0.05, 0.1], (0.05, 0.1)] ) def test_results_for_alpha_as_float_and_arraylike( - strategy: str, alpha: Any + estimator: BaseEstimator, strategy: str, alpha: Any ) -> None: """Test that output values do not depend on type of alpha.""" - mapie_reg = MapieRegressor(**STRATEGIES[strategy]) + mapie_reg = estimator(**STRATEGIES[strategy]) mapie_reg.fit(X, y) y_pred_float1, y_pis_float1 = mapie_reg.predict(X, alpha=alpha[0]) y_pred_float2, y_pis_float2 = mapie_reg.predict(X, alpha=alpha[1]) @@ -305,27 +386,33 @@ def test_results_for_alpha_as_float_and_arraylike( np.testing.assert_allclose(y_pis_float2[:, :, 0], y_pis_array[:, :, 1]) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_results_for_ordered_alpha(strategy: str) -> None: +def test_results_for_ordered_alpha( + estimator: BaseEstimator, strategy: str +) -> None: """ Test that prediction intervals lower (upper) bounds give consistent results for ordered alphas. """ - mapie = MapieRegressor(**STRATEGIES[strategy]) - mapie.fit(X, y) - y_pred, y_pis = mapie.predict(X, alpha=[0.05, 0.1]) + mapie_reg = estimator(**STRATEGIES[strategy]) + mapie_reg.fit(X, y) + _, y_pis = mapie_reg.predict(X, alpha=[0.05, 0.1]) assert (y_pis[:, 0, 0] <= y_pis[:, 0, 1]).all() assert (y_pis[:, 1, 0] >= y_pis[:, 1, 1]).all() +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_results_single_and_multi_jobs(strategy: str) -> None: +def test_results_single_and_multi_jobs( + estimator: BaseEstimator, strategy: str +) -> None: """ Test that MapieRegressor gives equal predictions regardless of number of parallel jobs. """ - mapie_single = MapieRegressor(n_jobs=1, **STRATEGIES[strategy]) - mapie_multi = MapieRegressor(n_jobs=-1, **STRATEGIES[strategy]) + mapie_single = estimator(n_jobs=1, **STRATEGIES[strategy]) + mapie_multi = estimator(n_jobs=-1, **STRATEGIES[strategy]) mapie_single.fit(X_toy, y_toy) mapie_multi.fit(X_toy, y_toy) y_pred_single, y_pis_single = mapie_single.predict(X_toy, alpha=0.2) @@ -334,16 +421,19 @@ def test_results_single_and_multi_jobs(strategy: str) -> None: np.testing.assert_allclose(y_pis_single, y_pis_multi) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_results_with_constant_sample_weights(strategy: str) -> None: +def test_results_with_constant_sample_weights( + estimator: BaseEstimator, strategy: str +) -> None: """ Test predictions when sample weights are None or constant with different values. """ n_samples = len(X) - mapie0 = MapieRegressor(**STRATEGIES[strategy]) - mapie1 = MapieRegressor(**STRATEGIES[strategy]) - mapie2 = MapieRegressor(**STRATEGIES[strategy]) + mapie0 = estimator(**STRATEGIES[strategy]) + mapie1 = estimator(**STRATEGIES[strategy]) + mapie2 = estimator(**STRATEGIES[strategy]) mapie0.fit(X, y, sample_weight=None) mapie1.fit(X, y, sample_weight=np.ones(shape=n_samples)) mapie2.fit(X, y, sample_weight=np.ones(shape=n_samples) * 5) @@ -356,24 +446,30 @@ def test_results_with_constant_sample_weights(strategy: str) -> None: np.testing.assert_allclose(y_pis1, y_pis2) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_prediction_between_low_up(strategy: str) -> None: +def test_prediction_between_low_up( + estimator: BaseEstimator, strategy: str +) -> None: """Test that prediction lies between low and up prediction intervals.""" - mapie = MapieRegressor(**STRATEGIES[strategy]) + mapie = estimator(**STRATEGIES[strategy]) mapie.fit(X, y) y_pred, y_pis = mapie.predict(X, alpha=0.1) assert (y_pred >= y_pis[:, 0, 0]).all() assert (y_pred <= y_pis[:, 1, 0]).all() +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("method", ["plus", "minmax"]) @pytest.mark.parametrize("agg_function", ["mean", "median"]) -def test_prediction_agg_function(method: str, agg_function: str) -> None: +def test_prediction_agg_function( + estimator: BaseEstimator, method: str, agg_function: str +) -> None: """ Test that predictions differ when ensemble is True/False, but not prediction intervals. """ - mapie = MapieRegressor(method=method, cv=2, agg_function=agg_function) + mapie = estimator(method=method, cv=2, agg_function=agg_function) mapie.fit(X, y) y_pred_1, y_pis_1 = mapie.predict(X, ensemble=True, alpha=0.1) y_pred_2, y_pis_2 = mapie.predict(X, ensemble=False, alpha=0.1) @@ -383,44 +479,37 @@ def test_prediction_agg_function(method: str, agg_function: str) -> None: np.testing.assert_allclose(y_pred_1, y_pred_2) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_linear_data_confidence_interval(strategy: str) -> None: - """ - Test that MapieRegressor applied on a linear regression model - fitted on a linear curve results in null uncertainty. - """ - mapie = MapieRegressor(**STRATEGIES[strategy]) - mapie.fit(X_toy, y_toy) - y_pred, y_pis = mapie.predict(X_toy, alpha=0.2) - np.testing.assert_allclose(y_pis[:, 0, 0], y_pis[:, 1, 0]) - np.testing.assert_allclose(y_pred, y_pis[:, 0, 0]) - - -@pytest.mark.parametrize("strategy", [*STRATEGIES]) -def test_linear_regression_results(strategy: str) -> None: +def test_linear_regression_results( + estimator: BaseEstimator, strategy: str +) -> None: """ Test expected prediction intervals for a multivariate linear regression problem with fixed random state. """ - mapie = MapieRegressor(**STRATEGIES[strategy]) + mapie = estimator(**STRATEGIES[strategy]) mapie.fit(X, y) _, y_pis = mapie.predict(X, alpha=0.05) y_pred_low, y_pred_up = y_pis[:, 0, 0], y_pis[:, 1, 0] width_mean = (y_pred_up - y_pred_low).mean() coverage = regression_coverage_score(y, y_pred_low, y_pred_up) - np.testing.assert_allclose(width_mean, WIDTHS[strategy], rtol=1e-2) - np.testing.assert_allclose(coverage, COVERAGES[strategy], rtol=1e-2) + np.testing.assert_allclose( + width_mean, WIDTHS[estimator.__name__][strategy], rtol=1e-2 + ) + np.testing.assert_allclose( + coverage, COVERAGES[estimator.__name__][strategy], rtol=1e-2 + ) -def test_results_prefit_ignore_method() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_results_prefit_ignore_method(estimator: BaseEstimator) -> None: """Test that method is ignored when ``cv="prefit"``.""" - estimator = LinearRegression().fit(X, y) + model = PREFIT_ESTIMATORS[estimator.__name__] all_y_pis: List[NDArray] = [] for method in METHODS: - mapie_reg = MapieRegressor( - estimator=estimator, cv="prefit", method=method - ) + mapie_reg = estimator(estimator=model, cv="prefit", method=method) mapie_reg.fit(X, y) _, y_pis = mapie_reg.predict(X, alpha=0.1) all_y_pis.append(y_pis) @@ -428,88 +517,46 @@ def test_results_prefit_ignore_method() -> None: np.testing.assert_allclose(y_pis1, y_pis2) -def test_results_prefit_naive() -> None: - """ - Test that prefit, fit and predict on the same dataset - is equivalent to the "naive" method. - """ - estimator = LinearRegression().fit(X, y) - mapie_reg = MapieRegressor(estimator=estimator, cv="prefit") - mapie_reg.fit(X, y) - _, y_pis = mapie_reg.predict(X, alpha=0.05) - width_mean = (y_pis[:, 1, 0] - y_pis[:, 0, 0]).mean() - coverage = regression_coverage_score(y, y_pis[:, 0, 0], y_pis[:, 1, 0]) - np.testing.assert_allclose(width_mean, WIDTHS["naive"], rtol=1e-2) - np.testing.assert_allclose(coverage, COVERAGES["naive"], rtol=1e-2) - - -def test_results_prefit() -> None: - """Test prefit results on a standard train/validation/test split.""" - X_train_val, X_test, y_train_val, y_test = train_test_split( - X, y, test_size=1 / 10, random_state=1 - ) - X_train, X_val, y_train, y_val = train_test_split( - X_train_val, y_train_val, test_size=1 / 9, random_state=1 - ) - estimator = LinearRegression().fit(X_train, y_train) - mapie_reg = MapieRegressor(estimator=estimator, cv="prefit") - mapie_reg.fit(X_val, y_val) - _, y_pis = mapie_reg.predict(X_test, alpha=0.05) - width_mean = (y_pis[:, 1, 0] - y_pis[:, 0, 0]).mean() - coverage = regression_coverage_score( - y_test, y_pis[:, 0, 0], y_pis[:, 1, 0] - ) - np.testing.assert_allclose(width_mean, WIDTHS["prefit"], rtol=1e-2) - np.testing.assert_allclose(coverage, COVERAGES["prefit"], rtol=1e-2) - - -def test_not_enough_resamplings() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_not_enough_resamplings(estimator: BaseEstimator) -> None: """ Test that a warning is raised if at least one conformity score is nan. """ with pytest.warns(UserWarning, match=r"WARNING: at least one point of*"): - mapie_reg = MapieRegressor( + mapie_reg = estimator( cv=Subsample(n_resamplings=1), agg_function="mean" ) mapie_reg.fit(X, y) -def test_no_agg_fx_specified_with_subsample() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_no_agg_fx_specified_with_subsample(estimator: BaseEstimator) -> None: """ Test that a warning is raised if at least one conformity score is nan. """ with pytest.raises( ValueError, match=r"You need to specify an aggregation*" ): - mapie_reg = MapieRegressor( + mapie_reg = estimator( cv=Subsample(n_resamplings=1), agg_function=None ) mapie_reg.fit(X, y) -def test_invalid_aggregate_all() -> None: - """ - Test that wrong aggregation in MAPIE raise errors. - """ - with pytest.raises( - ValueError, - match=r".*Aggregation function called but not defined.*", - ): - aggregate_all(None, X) - - -def test_aggregate_with_mask_with_prefit() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_aggregate_with_mask_with_prefit(estimator) -> None: """ Test ``_aggregate_with_mask`` in case ``cv`` is ``"prefit"``. """ - mapie_reg = MapieRegressor(cv="prefit") + k = None + mapie_reg = estimator(cv="prefit") with pytest.raises( ValueError, match=r".*There should not be aggregation of predictions if cv is*", ): mapie_reg._aggregate_with_mask(k, k) - mapie_reg = MapieRegressor(agg_function="nonsense") + mapie_reg = estimator(agg_function="nonsense") with pytest.raises( ValueError, match=r".*The value of self.agg_function is not correct*", @@ -517,12 +564,12 @@ def test_aggregate_with_mask_with_prefit() -> None: mapie_reg._aggregate_with_mask(k, k) -def test_pred_loof_isnan() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_pred_oof_isnan(estimator: BaseEstimator) -> None: """Test that if validation set is empty then prediction is empty.""" - mapie_reg = MapieRegressor() - y_pred: ArrayLike + mapie_reg = estimator() _, y_pred, _ = mapie_reg._fit_and_predict_oof_model( - estimator=LinearRegression(), + estimator=BASE_ESTIMATORS[estimator.__name__](), X=X_toy, y=y_toy, train_index=[0, 1, 2, 3, 4], @@ -531,16 +578,22 @@ def test_pred_loof_isnan() -> None: assert len(y_pred) == 0 -def test_pipeline_compatibility() -> None: +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) +def test_pipeline_compatibility(estimator: BaseEstimator) -> None: """Check that MAPIE works on pipeline based on pandas dataframes""" + X_toy, y_toy = make_regression( + n_samples=100, n_features=1, noise=1.0, random_state=1 + ) X = pd.DataFrame( { - "x_cat": ["A", "A", "B", "A", "A", "B"], - "x_num": [0, 1, 1, 4, np.nan, 5], - "y": [5, 7, 3, 9, 10, 8], + "x_cat": np.random.choice( + ["A", "B"], X_toy.shape[0], replace=True + ), + "x_num": X_toy.flatten(), + "y": y_toy, } ) - y = pd.Series([5, 7, 3, 9, 10, 8]) + y = pd.Series(y_toy) numeric_preprocessor = Pipeline( [ ("imputer", SimpleImputer(strategy="mean")), @@ -555,21 +608,22 @@ def test_pipeline_compatibility() -> None: ("num", numeric_preprocessor, ["x_num"]), ] ) - pipe = make_pipeline(preprocessor, LinearRegression()) - mapie = MapieRegressor(pipe) + pipe = make_pipeline(preprocessor, BASE_ESTIMATORS[estimator.__name__]()) + mapie = estimator(pipe) mapie.fit(X, y) mapie.predict(X) +@pytest.mark.parametrize("estimator", MAPIE_ESTIMATORS) @pytest.mark.parametrize("strategy", [*STRATEGIES]) @pytest.mark.parametrize( "conformity_score", [AbsoluteConformityScore(), GammaConformityScore()] ) def test_gammaconformityscore( - strategy: str, conformity_score: ConformityScore + estimator: BaseEstimator, strategy: str, conformity_score: ConformityScore ) -> None: """Test that GammaConformityScore with MAPIE raises no error.""" - mapie_reg = MapieRegressor( + mapie_reg = estimator( conformity_score=conformity_score, **STRATEGIES[strategy] ) diff --git a/mapie/tests/test_utils.py b/mapie/tests/test_utils.py index d7af1ecfd..30a48d4cb 100644 --- a/mapie/tests/test_utils.py +++ b/mapie/tests/test_utils.py @@ -324,33 +324,6 @@ def test_compute_quantiles_2D_and_3D(alphas: NDArray): assert (quantiles1 == quantiles2).all() -@pytest.mark.parametrize("estimator", [-1, 3, 0.2]) -def test_quantile_prefit_non_iterable(estimator: Any) -> None: - """ - Test that there is a list of estimators provided when cv='prefit' - is called for MapieQuantileRegressor. - """ - with pytest.raises( - ValueError, - match=r".*Estimator for prefit must be an iterable object.*", - ): - mapie_reg = MapieQuantileRegressor(estimator=estimator, cv="prefit") - mapie_reg.fit([1, 2, 3], [4, 5, 6]) - - -# def test_calib_set_no_Xy_but_sample_weight() -> None: -# """Test warning message if sample weight provided but no X y in calib.""" -# X = np.array([4, 5, 6]) -# y = np.array([4, 3, 2]) -# sample_weight = np.array([4, 4, 4]) -# sample_weight_calib = np.array([4, 3, 4]) -# with pytest.warns(UserWarning, match=r"WARNING: sample weight*"): -# check_calib_set( -# X=X, y=y, sample_weight=sample_weight, -# sample_weight_calib=sample_weight_calib -# ) - - @pytest.mark.parametrize("strategy", ["quantile", "uniform", "array split"]) def test_binning_group_strategies(strategy: str) -> None: """Test that different strategies have the correct outputs.""" diff --git a/mapie/utils.py b/mapie/utils.py index bfd946f50..67150ad2e 100644 --- a/mapie/utils.py +++ b/mapie/utils.py @@ -321,7 +321,7 @@ def check_n_features_in( def check_alpha_and_n_samples( - alphas: Union[Iterable[float], float], + alphas: Iterable[float], n: int, ) -> None: """ @@ -354,8 +354,6 @@ def check_alpha_and_n_samples( Number of samples of the score is too low, 1/alpha (or 1/(1 - alpha)) must be lower than the number of samples. """ - if isinstance(alphas, float): - alphas = np.array([alphas]) for alpha in alphas: if n < 1 / alpha or n < 1 / (1 - alpha): raise ValueError(