From 85f37e9c56e94474fc31fb2f19b0ee51fb54987f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Henrik=20Bostr=C3=B6m?= Date: Sat, 21 Sep 2024 17:33:49 +0200 Subject: [PATCH] 0.7.1 --- src/crepes/base.py | 1514 ++++++++++++++++++++++++------------------ src/crepes/extras.py | 155 +++-- 2 files changed, 959 insertions(+), 710 deletions(-) diff --git a/src/crepes/base.py b/src/crepes/base.py index bf10392..98766f9 100644 --- a/src/crepes/base.py +++ b/src/crepes/base.py @@ -14,7 +14,7 @@ """ -__version__ = "0.7.0" +__version__ = "0.7.1" import numpy as np import pandas as pd @@ -39,6 +39,7 @@ def __init__(self): self.time_evaluate = None self.alphas = None self.normalized = None + self.seed = None class ConformalClassifier(ConformalPredictor): """ @@ -53,7 +54,7 @@ def __repr__(self): else: return f"ConformalClassifier(fitted={self.fitted})" - def fit(self, alphas, bins=None): + def fit(self, alphas, bins=None, seed=None): """ Fit conformal classifier. @@ -63,6 +64,8 @@ def fit(self, alphas, bins=None): non-conformity scores bins : array-like of shape (n_samples,), default=None Mondrian categories + seed : int, default=None + set random seed Returns ------- @@ -90,6 +93,12 @@ def fit(self, alphas, bins=None): cc_mond = ConformalClassifier() cc_mond.fit(alphas_cal, bins=bins_cal) + + Note + ---- + By providing a random seed, e.g., ``seed=123``, calls to the methods + ``predict_p``, ``predict_set`` and ``evaluate`` of the + :class:`.ConformalClassifier` object will be deterministic. """ tic = time.time() if bins is None: @@ -100,14 +109,16 @@ def fit(self, alphas, bins=None): bin_values = np.unique(bins) self.alphas = (bin_values, [np.sort(alphas[bins==b])[::-1] for b in bin_values]) + self.seed = seed self.fitted = True toc = time.time() self.time_fit = toc-tic return self - def predict_p(self, alphas, bins=None, confidence=0.95): + def predict_p(self, alphas, bins=None, confidence=0.95, smoothing=True, + seed=None): """ - Obtain (smoothed) p-values from conformal classifier. + Obtain (smoothed or non-smoothed) p-values from conformal classifier. Parameters ---------- @@ -117,6 +128,10 @@ def predict_p(self, alphas, bins=None, confidence=0.95): Mondrian categories confidence : float in range (0,1), default=0.95 confidence level + smoothing : bool, default=True + use smoothed p-values + seed : int, default=None + set random seed Returns ------- @@ -135,33 +150,60 @@ def predict_p(self, alphas, bins=None, confidence=0.95): Assuming that ``bins_test`` is a vector with Mondrian categories (bin labels) for the test set and ``cc_mond`` a fitted Mondrian conformal - classifier, then the following provides p-values for the test set: + classifier, then the following provides (smoothed) p-values for the + test set: .. code-block:: python p_values = cc_mond.predict_p(alphas_test, bins=bins_test) + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``fit``. """ tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) if not self.mondrian: - p_values = np.array( - [[(np.sum(self.alphas > alpha) + np.random.rand()*( - np.sum(self.alphas == alpha)+1))/(len(self.alphas)+1) - for alpha in alpha_row] for alpha_row in alphas]) - else: + if smoothing: + p_values = np.array( + [[(np.sum(self.alphas > alpha) + np.random.rand()*( + np.sum(self.alphas == alpha)+1))/(len(self.alphas)+1) + for alpha in alpha_row] for alpha_row in alphas]) + else: + p_values = np.array( + [[(np.sum(self.alphas >= alpha) + 1)/(len(self.alphas)+1) + for alpha in alpha_row] for alpha_row in alphas]) + else: bin_values, bin_alphas = self.alphas bin_indexes = np.array([np.argwhere(bin_values == bins[i])[0][0] for i in range(len(bins))]) - p_values = np.array([[(np.sum(bin_alphas[bin_indexes[i]] > alpha) \ - + np.random.rand()*(np.sum(bin_alphas[ - bin_indexes[i]] == alpha)+1))/( - len(bin_alphas[bin_indexes[i]])+1) - for alpha in alphas[i]] - for i in range(len(alphas))]) + if smoothing: + p_values = np.array([ + [(np.sum(bin_alphas[bin_indexes[i]] > alpha) \ + + np.random.rand()*(np.sum(bin_alphas[ + bin_indexes[i]] == alpha)+1))/( + len(bin_alphas[bin_indexes[i]])+1) + for alpha in alphas[i]] + for i in range(len(alphas))]) + else: + p_values = np.array([ + [(np.sum(bin_alphas[bin_indexes[i]] >= alpha) + 1)/( + len(bin_alphas[bin_indexes[i]])+1) + for alpha in alphas[i]] + for i in range(len(alphas))]) + if seed is not None: + np.random.set_state(random_state) toc = time.time() self.time_predict = toc-tic return p_values - def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=False): + def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=True, + seed=None): """ Obtain prediction sets using conformal classifier. @@ -173,8 +215,10 @@ def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=False): Mondrian categories confidence : float in range (0,1), default=0.95 confidence level - smoothing : bool, default=False + smoothing : bool, default=True use smoothed p-values + seed : int, default=None + set random seed Returns ------- @@ -205,13 +249,23 @@ def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=False): Note ---- - Using smoothed p-values substantially increases computation time and - hardly has any effect on the predictions sets, except for when having - small calibration sets. + The use of smoothed p-values increases computation time and typically + has a minor effect on the predictions sets, except for small calibration + sets. + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``fit``. """ tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) if smoothing: - p_values = self.predict_p(alphas, bins) + p_values = self.predict_p(alphas, bins, smoothing=True) prediction_sets = (p_values >= 1-confidence).astype(int) elif bins is None: alpha_index = int((1-confidence)*(len(self.alphas)+1))-1 @@ -241,12 +295,14 @@ def predict_set(self, alphas, bins=None, confidence=0.95, smoothing=False): " too small for the chosen confidence level; " \ "all labels are included in the corresponding" \ "prediction sets") + if seed is not None: + np.random.set_state(random_state) toc = time.time() self.time_predict = toc-tic return prediction_sets def evaluate(self, alphas, classes, y, bins=None, confidence=0.95, - smoothing=False, metrics=None): + smoothing=True, metrics=None, seed=None): """ Evaluate conformal classifier. @@ -262,11 +318,13 @@ class names Mondrian categories confidence : float in range (0,1), default=0.95 confidence level - smoothing : bool, default=False + smoothing : bool, default=True use smoothed p-values metrics : a string or a list of strings, default = list of all metrics, i.e., ["error", "avg_c", "one_c", "empty", "time_fit", "time_evaluate"] + seed : int, default=None + set random seed Returns ------- @@ -294,16 +352,27 @@ class names Note ---- - Using smoothed p-values substantially increases computation time and - hardly has any effect on the results, except for when having small - calibration sets. + The use of smoothed p-values increases computation time and typically + has a minor effect on the results, except for small calibration sets. + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``fit``. """ if metrics is None: metrics = ["error", "avg_c", "one_c", "empty", "time_fit", "time_evaluate"] tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) prediction_sets = self.predict_set(alphas, bins, confidence, smoothing) test_results = get_test_results(prediction_sets, classes, y, metrics) + if seed is not None: + np.random.set_state(random_state) toc = time.time() self.time_evaluate = toc-tic if "time_fit" in metrics: @@ -633,7 +702,7 @@ def __repr__(self): else: return f"ConformalPredictiveSystem(fitted={self.fitted})" - def fit(self, residuals, sigmas=None, bins=None): + def fit(self, residuals, sigmas=None, bins=None, seed=None): """ Fit conformal predictive system. @@ -645,6 +714,8 @@ def fit(self, residuals, sigmas=None, bins=None): difficulty estimates bins : array-like of shape (n_values,), default=None Mondrian categories + seed : int, default=None + set random seed Returns ------- @@ -693,8 +764,13 @@ def fit(self, residuals, sigmas=None, bins=None): cps_norm_mond = ConformalPredictiveSystem() cps_norm_mond.fit(residuals_cal, sigmas=sigmas_cal, bins=bins_cal) - """ + Note + ---- + By providing a random seed, e.g., ``seed=123``, calls to the methods + ``predict`` and ``evaluate`` of the :class:`.ConformalPredictiveSystem` + object will be deterministic. + """ tic = time.time() if bins is None: self.mondrian = False @@ -716,6 +792,7 @@ def fit(self, residuals, sigmas=None, bins=None): self.alphas = (bin_values, [np.sort( residuals[bins==b]/sigmas[bins==b]) for b in bin_values]) self.fitted = True + self.seed = seed toc = time.time() self.time_fit = toc-tic return self @@ -723,7 +800,7 @@ def fit(self, residuals, sigmas=None, bins=None): def predict(self, y_hat, sigmas=None, bins=None, y=None, lower_percentiles=None, higher_percentiles=None, y_min=-np.inf, y_max=np.inf, return_cpds=False, - cpds_by_bins=False): + cpds_by_bins=False, seed=None): """ Predict using conformal predictive system. @@ -755,6 +832,8 @@ def predict(self, y_hat, sigmas=None, bins=None, cpds_by_bins : Boolean, default=False specifies whether the output cpds should be grouped by bin or not; only applicable when bins is not None and return_cpds = True + seed : int, default=None + set random seed Returns ------- @@ -861,9 +940,18 @@ def predict(self, y_hat, sigmas=None, bins=None, ---- Setting ``cpds_by_bins=True`` has an effect only for Mondrian conformal predictive systems. - """ + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``fit``. + """ tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) if self.mondrian: bin_values, bin_alphas = self.alphas bin_indexes = [np.argwhere(bins == b).T[0] for b in bin_values] @@ -1170,9 +1258,12 @@ def predict(self, y_hat, sigmas=None, bins=None, cpds_out[bin_indexes[b]] = [cpds[b][i] for i in range(len(cpds[b]))] return cpds_out + if seed is not None: + np.random.set_state(random_state) def evaluate(self, y_hat, y, sigmas=None, bins=None, - confidence=0.95, y_min=-np.inf, y_max=np.inf, metrics=None): + confidence=0.95, y_min=-np.inf, y_max=np.inf, + metrics=None, seed=None): """ Evaluate conformal predictive system. @@ -1195,6 +1286,8 @@ def evaluate(self, y_hat, y, sigmas=None, bins=None, metrics : a string or a list of strings, default=list of all metrics; ["error", "eff_mean","eff_med", "CRPS", "time_fit", "time_evaluate"] + seed : int, default=None + set random seed Returns ------- @@ -1226,12 +1319,22 @@ def evaluate(self, y_hat, y, sigmas=None, bins=None, number of calibration and test objects, unless a Mondrian approach is employed; for the latter, this number is reduced by increasing the number of bins. + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``fit``. """ + tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) if isinstance(y, pd.Series): y = y.values if isinstance(y_hat, pd.Series): y_hat = y_hat.values - tic = time.time() test_results = {} lower_percentile = (1-confidence)/2*100 higher_percentile = (confidence+(1-confidence)/2)*100 @@ -1287,7 +1390,9 @@ def evaluate(self, y_hat, y, sigmas=None, bins=None, if "time_fit" in metrics: test_results["time_fit"] = self.time_fit toc = time.time() - self.time_evaluate = toc-tic + if seed is not None: + np.random.set_state(random_state) + self.time_evaluate = toc-tic if "time_evaluate" in metrics: test_results["time_evaluate"] = self.time_evaluate return test_results @@ -1367,32 +1472,25 @@ def get_crps(cpd_index, lower_errors, higher_errors, widths, sigma, cpd, y): higher_errors[cpd_index]*(cpd[cpd_index+1]-y)*sigma return score -class WrapRegressor(): +class WrapClassifier(): """ - A learner wrapped with a :class:`.ConformalRegressor` - or :class:`.ConformalPredictiveSystem`. + A learner wrapped with a :class:`.ConformalClassifier`. """ def __init__(self, learner): - self.cr = None - self.cps = None + self.cc = None + self.nc = None self.calibrated = False self.learner = learner - self.de = None - self.mc = None - + self.seed = None + def __repr__(self): if self.calibrated: - if self.cr is not None: - return (f"WrapRegressor(learner={self.learner}, " - f"calibrated={self.calibrated}, " - f"predictor={self.cr})") - else: - return (f"WrapRegressor(learner={self.learner}, " - f"calibrated={self.calibrated}, " - f"predictor={self.cps})") + return (f"WrapClassifier(learner={self.learner}, " + f"calibrated={self.calibrated}, " + f"predictor={self.cc})") else: - return f"WrapRegressor(learner={self.learner}, calibrated={self.calibrated})" + return f"WrapClassifier(learner={self.learner}, calibrated={self.calibrated})" def fit(self, X, y, **kwargs): """ @@ -1420,10 +1518,10 @@ def fit(self, X, y, **kwargs): .. code-block:: python - from sklearn.ensemble import RandomForestRegressor - from crepes import WrapRegressor + from sklearn.ensemble import RandomForestClassifier + from crepes import WrapClassifier - rf = WrapRegressor(RandomForestRegressor()) + rf = Wrap(RandomForestClassifier()) rf.fit(X_train, y_train) Note @@ -1448,7 +1546,7 @@ def predict(self, X): Parameters ---------- X : array-like of shape (n_samples, n_features), - set of objects + set of objects Returns ------- @@ -1458,9 +1556,10 @@ def predict(self, X): Examples -------- - Assuming ``w`` is a :class:`.WrapRegressor` object for which the wrapped - learner ``w.learner`` has been fitted, (point) predictions of the - learner can be obtained for a set of test objects ``X_test`` by: + Assuming ``w`` is a :class:`.WrapClassifier` object for which the + wrapped learner ``w.learner`` has been fitted, (point) + predictions of the learner can be obtained for a set + of test objects ``X_test`` by: .. code-block:: python @@ -1474,10 +1573,44 @@ def predict(self, X): """ return self.learner.predict(X) - def calibrate(self, X, y, de=None, mc=None, oob=False, cps=False): + def predict_proba(self, X): """ - Fit a :class:`.ConformalRegressor` or - :class:`.ConformalPredictiveSystem` using the wrapped learner. + Predict with learner. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features), + set of objects + + Returns + ------- + y : array-like of shape (n_samples, n_classes), + predicted probabilities using the ``predict_proba`` + method of the ``learner`` object. + + Examples + -------- + Assuming ``w`` is a :class:`.WrapClassifier` object for which the + wrapped learner ``w.learner`` has been fitted, predicted + probabilities of the learner can be obtained for a set + of test objects ``X_test`` by: + + .. code-block:: python + + probabilities = w.predict_proba(X_test) + + The above is equivalent to: + + .. code-block:: python + + probabilities = w.learner.predict_proba(X_test) + """ + return self.learner.predict_proba(X) + + def calibrate(self, X, y, oob=False, class_cond=False, nc=hinge, mc=None, + seed=None): + """ + Fit a :class:`.ConformalClassifier` using the wrapped learner. Parameters ---------- @@ -1485,55 +1618,45 @@ def calibrate(self, X, y, de=None, mc=None, oob=False, cps=False): set of objects y : array-like of shape (n_samples,), target values - de: :class:`crepes.extras.DifficultyEstimator`, default=None - object used for computing difficulty estimates - mc: function or :class:`crepes.extras.MondrianCategorizer`, default=None - function or :class:`crepes.extras.MondrianCategorizer` for computing Mondrian categories oob : bool, default=False use out-of-bag estimation - cps : bool, default=False - if cps=False, the method fits a :class:`.ConformalRegressor` - and otherwise, a :class:`.ConformalPredictiveSystem` + class_cond : bool, default=False + if class_cond=True, the method fits a Mondrian + :class:`.ConformalClassifier` using the class + labels as categories + nc : function, default = :func:`crepes.extras.hinge` + function to compute non-conformity scores + mc: function or :class:`crepes.extras.MondrianCategorizer`, default=None + function or :class:`crepes.extras.MondrianCategorizer` for computing Mondrian + categories + seed : int, default=None + set random seed Returns ------- self : object - The :class:`.WrapRegressor` object is updated with a fitted - :class:`.ConformalRegressor` or :class:`.ConformalPredictiveSystem` + Wrap object updated with a fitted :class:`.ConformalClassifier` Examples -------- Assuming ``X_cal`` and ``y_cal`` to be an array and vector, respectively, with objects and labels for the calibration set, - and ``w`` is a :class:`.WrapRegressor` object for which the learner - has been fitted, a standard conformal regressor is formed by: + and ``w`` is a :class:`.WrapClassifier` object for which the learner + has been fitted, a standard conformal classifier can be formed by: .. code-block:: python w.calibrate(X_cal, y_cal) - Assuming that ``de`` is a fitted difficulty estimator, - a normalized conformal regressor is obtained by: - - .. code-block:: python - - w.calibrate(X_cal, y_cal, de=de) - - Assuming that ``get_categories`` is a function that returns categories - (bin labels), a Mondrian conformal regressor is obtained by: + Assuming that ``get_categories`` is a function that returns a vector of + Mondrian categories (bin labels), a Mondrian conformal classifier can + be generated by: .. code-block:: python w.calibrate(X_cal, y_cal, mc=get_categories) - A normalized Mondrian conformal regressor is generated in the - following way: - - .. code-block:: python - - w.calibrate(X_cal, y_cal, de=de, mc=get_categories) - - By providing the option ``oob=True``, the conformal regressor + By providing the option ``oob=True``, the conformal classifier will be calibrating using out-of-bag predictions, allowing the full set of training objects (``X_train``) and labels (``y_train``) to be used, e.g., @@ -1542,357 +1665,597 @@ def calibrate(self, X, y, de=None, mc=None, oob=False, cps=False): w.calibrate(X_train, y_train, oob=True) - By providing the option ``cps=True``, each of the above calls will instead - generate a :class:`.ConformalPredictiveSystem`, e.g., + By providing the option ``class_cond=True``, a Mondrian conformal classifier + will be formed using the class labels as categories, e.g., .. code-block:: python - w.calibrate(X_cal, y_cal, de=de, cps=True) + w.calibrate(X_cal, y_cal, class_cond=True) + + Note + ---- + Any Mondrian categorizer specified by the ``mc`` argument will be + ignored by :meth:`.calibrate`, if ``class_cond=True``, as the latter + implies that Mondrian categories are formed using the labels in ``y``. + + Note + ---- + By providing a random seed, e.g., ``seed=123``, the call to ``calibrate`` + as well as calls to the methods ``predict_set``, ``predict_p`` and + ``evaluate`` of the :class:`.WrapClassifier` object will be + deterministic. Note ---- Enabling out-of-bag calibration, i.e., setting ``oob=True``, requires - that the wrapped learner has an attribute ``oob_prediction_``, which - e.g., is the case for a ``sklearn.ensemble.RandomForestRegressor``, if - enabled when created, e.g., ``RandomForestRegressor(oob_score=True)`` + that the wrapped learner has an attribute ``oob_decision_function_``, + which e.g., is the case for a ``sklearn.ensemble.RandomForestClassifier``, + if enabled when created, e.g., ``RandomForestClassifier(oob_score=True)`` Note ---- - The use of out-of-bag calibration, as enabled by ``oob=True``, - does not come with the theoretical validity guarantees of the regular - (inductive) conformal regressors and predictive systems, due to that - calibration and test instances are not handled in exactly the same way. + The use of out-of-bag calibration, as enabled by ``oob=True``, does not + come with the theoretical validity guarantees of the regular (inductive) + conformal classifiers, due to that calibration and test instances are not + handled in exactly the same way. """ + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + self.seed = seed if isinstance(y, pd.Series): y = y.values + self.cc = ConformalClassifier() + self.nc = nc + self.mc = mc + self.class_cond = class_cond if oob: - residuals = y - self.learner.oob_prediction_ - else: - residuals = y - self.predict(X) - if de is None: - sigmas = None - else: - sigmas = de.apply(X) - self.de = de - if mc is None: - bins = None - elif isinstance(mc, MondrianCategorizer): - bins = mc.apply(X) + alphas = nc(self.learner.oob_decision_function_, self.learner.classes_, y) else: - bins = mc(X) - self.mc = mc - if not cps: - self.cr = ConformalRegressor() - self.cr.fit(residuals, sigmas=sigmas, bins=bins) - self.cps = None + alphas = nc(self.learner.predict_proba(X), self.learner.classes_, y) + if class_cond: + self.cc.fit(alphas, bins=y) else: - self.cps = ConformalPredictiveSystem() - self.cps.fit(residuals, sigmas=sigmas, bins=bins) - self.cr = None + if isinstance(mc, MondrianCategorizer): + bins = mc.apply(X) + self.cc.fit(alphas, bins=bins) + elif mc is not None: + bins = mc(X) + self.cc.fit(alphas, bins=bins) + else: + self.cc.fit(alphas) self.calibrated = True + if seed is not None: + np.random.set_state(random_state) return self - def predict_int(self, X, confidence=0.95, y_min=-np.inf, y_max=np.inf): + def predict_p(self, X, smoothing=True, seed=None): """ - Predict interval using fitted :class:`.ConformalRegressor` or - :class:`.ConformalPredictiveSystem`. + Obtain (smoothed or non-smoothed) p-values using conformal classifier. Parameters ---------- X : array-like of shape (n_samples, n_features), set of objects - confidence : float in range (0,1), default=0.95 - confidence level - y_min : float or int, default=-numpy.inf - minimum value to include in prediction intervals - y_max : float or int, default=numpy.inf - maximum value to include in prediction intervals + smoothing : bool, default=True + use smoothed p-values + seed : int, default=None + set random seed Returns ------- - intervals : ndarray of shape (n_samples, 2) - prediction intervals + p-values : ndarray of shape (n_samples, n_classes) + p-values Examples -------- Assuming that ``X_test`` is a set of test objects and ``w`` is a - :class:`.WrapRegressor` object that has been calibrated, i.e., - :meth:`.calibrate` has been applied, prediction intervals at the - 99% confidence level can be obtained by: - - .. code-block:: python - - intervals = w.predict_int(X_test, confidence=0.99) - - The following provides prediction intervals at the default confidence - level (95%), where the intervals are lower-bounded by 0: + :class:`.WrapClassifier` object that has been calibrated, i.e., + :meth:`.calibrate` has been applied, the (smoothed) p-values for + the test objects are obtained by: .. code-block:: python - intervals = w.predict_int(X_test, y_min=0) + p_values = w.predict_p(X_test) Note ---- - In case the specified confidence level is too high in relation to the - size of the calibration set, a warning will be issued and the output - intervals will be of maximum size. + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. """ - if not self.calibrated: - raise RuntimeError(("predict_int requires that calibrate has been" - "called first")) - else: - y_hat = self.learner.predict(X) - if self.de is None: - sigmas = None - else: - sigmas = self.de.apply(X) - if self.mc is None: - bins = None - elif isinstance(self.mc, MondrianCategorizer): + tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + alphas = self.nc(self.learner.predict_proba(X)) + if self.class_cond: + p_values = np.array([ + self.cc.predict_p(alphas, + np.full(len(X), + self.learner.classes_[c]), + smoothing)[:, c] + for c in range(len(self.learner.classes_))]).T + else: + if isinstance(self.mc, MondrianCategorizer): bins = self.mc.apply(X) - else: + p_values = self.cc.predict_p(alphas, bins, smoothing) + elif self.mc is not None: bins = self.mc(X) - if self.cr is not None: - return self.cr.predict(y_hat, sigmas=sigmas, bins=bins, - confidence=confidence, - y_min=y_min, y_max=y_max) + p_values = self.cc.predict_p(alphas, bins, smoothing) else: - lower_percentile = (1-confidence)/2*100 - higher_percentile = (confidence+(1-confidence)/2)*100 - return self.cps.predict(y_hat, sigmas=sigmas, bins=bins, - lower_percentiles=lower_percentile, - higher_percentiles=higher_percentile, - y_min=y_min, y_max=y_max) + p_values = self.cc.predict_p(alphas, smoothing=smoothing) + if seed is not None: + np.random.set_state(random_state) + toc = time.time() + self.time_predict = toc-tic + return p_values - def predict_cps(self, X, y=None, lower_percentiles=None, - higher_percentiles=None, y_min=-np.inf, y_max=np.inf, - return_cpds=False, cpds_by_bins=False): + def predict_set(self, X, confidence=0.95, smoothing=True, seed=None): """ - Predict using :class:`.ConformalPredictiveSystem`. + Obtain prediction sets using conformal classifier. Parameters ---------- X : array-like of shape (n_samples, n_features), set of objects - y : float, int or array-like of shape (n_samples,), default=None - values for which p-values should be returned - lower_percentiles : array-like of shape (l_values,), default=None - percentiles for which a lower value will be output - in case a percentile lies between two values - (similar to `interpolation="lower"` in `numpy.percentile`) - higher_percentiles : array-like of shape (h_values,), default=None - percentiles for which a higher value will be output - in case a percentile lies between two values - (similar to `interpolation="higher"` in `numpy.percentile`) - y_min : float or int, default=-numpy.inf - The minimum value to include in prediction intervals. - y_max : float or int, default=numpy.inf - The maximum value to include in prediction intervals. - return_cpds : Boolean, default=False - specifies whether conformal predictive distributions (cpds) - should be output or not - cpds_by_bins : Boolean, default=False - specifies whether the output cpds should be grouped by bin or not; - only applicable when bins is not None and return_cpds = True + confidence : float in range (0,1), default=0.95 + confidence level + smoothing : bool, default=True + use smoothed p-values + seed : int, default=None + set random seed Returns ------- - results : ndarray of shape (n_samples, n_cols) or (n_samples,) - the shape is (n_samples, n_cols) if n_cols > 1 and otherwise - (n_samples,), where n_cols = p_values+l_values+h_values where - p_values = 1 if y is not None and 0 otherwise, l_values are the - number of lower percentiles, and h_values are the number of higher - percentiles. Only returned if n_cols > 0. - cpds : ndarray of (n_samples, c_values), ndarray of (n_samples,) - or list of ndarrays - conformal predictive distributions. Only returned if - return_cpds == True. For non-Mondrian conformal predictive systems, - the distributions are represented by a single array, where the - number of columns (c_values) is determined by the number of - residuals of the fitted conformal predictive system. For Mondrian - conformal predictive systems, the distributions are represented by - a vector of arrays, if cpds_by_bins = False, or a list of arrays, - with one element for each Mondrian category, if cpds_by_bins = True. + prediction sets : ndarray of shape (n_values, n_classes) + prediction sets, where the value 1 (0) indicates + that the class label is included (excluded), i.e., + the corresponding p-value is less than 1-confidence + + Examples + -------- + Assuming that ``X_test`` is a set of test objects and ``w`` is a + :class:`.WrapClassifier` object that has been calibrated, i.e., + :meth:`.calibrate` has been applied, the prediction sets for the + test objects at the 99% confidence level are obtained by: + + .. code-block:: python + + prediction_sets = w.predict_set(X_test, confidence=0.99) + + Note + ---- + The use of smoothed p-values increases computation time and typically + has a minor effect on the predictions sets, except for small calibration + sets. + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. + """ + tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + alphas = self.nc(self.learner.predict_proba(X)) + if self.class_cond: + prediction_set = np.array([ + self.cc.predict_set(alphas, + np.full(len(X), + self.learner.classes_[c]), + confidence, smoothing)[:, c] + for c in range(len(self.learner.classes_))]).T + else: + if isinstance(self.mc, MondrianCategorizer): + bins = self.mc.apply(X) + elif self.mc is not None: + bins = self.mc(X) + else: + bins = None + prediction_set = self.cc.predict_set(alphas, bins, confidence, + smoothing) + if seed is not None: + np.random.set_state(random_state) + toc = time.time() + self.time_predict = toc-tic + return prediction_set + + def evaluate(self, X, y, confidence=0.95, smoothing=True, + metrics=None, seed=None): + """ + Evaluate the conformal classifier. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features) + set of objects + y : array-like of shape (n_samples,) + correct target values + confidence : float in range (0,1), default=0.95 + confidence level + smoothing : bool, default=True + use smoothed p-values + metrics : a string or a list of strings, + default=list of all metrics, i.e., ["error", "avg_c", "one_c", + "empty", "time_fit", "time_evaluate"] + seed : int, default=None + set random seed + + Returns + ------- + results : dictionary with a key for each selected metric + estimated performance using the metrics, where "error" is the + fraction of prediction sets not containing the true class label, + "avg_c" is the average no. of predicted class labels, "one_c" is + the fraction of singleton prediction sets, "empty" is the fraction + of empty prediction sets, "time_fit" is the time taken to fit the + conformal classifier, and "time_evaluate" is the time taken for the + evaluation Examples -------- Assuming that ``X_test`` is a set of test objects, ``y_test`` is a - vector with true targets, ``w`` is a :class:`.WrapRegressor` object - calibrated with the option ``cps=True``, p-values for the true targets - can be obtained by: + vector with true targets, and ``w`` is a calibrated + :class:`.WrapClassifier` object, then the latter can be evaluated at + the 90% confidence level with respect to error, average prediction set + size and fraction of singleton predictions by: .. code-block:: python - p_values = w.predict_cps(X_test, y=y_test) + results = w.evaluate(X_test, y_test, confidence=0.9, + metrics=["error", "avg_c", "one_c"]) - P-values with respect to some specific value, e.g., 37, can be - obtained by: + Note + ---- + The reported result for ``time_fit`` only considers fitting the + conformal regressor or predictive system; not for fitting the + learner. + + Note + ---- + The use of smoothed p-values increases computation time and typically + has a minor effect on the results, except for small calibration sets. + + Note + ---- + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. + """ + if isinstance(y, pd.Series): + y = y.values + if not self.calibrated: + raise RuntimeError(("evaluate requires that calibrate has been" + "called first")) + else: + if metrics is None: + metrics = ["error", "avg_c", "one_c", "empty", "time_fit", + "time_evaluate"] + tic = time.time() + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + prediction_sets = self.predict_set(X, confidence, smoothing) + test_results = get_test_results(prediction_sets, + self.learner.classes_, y, metrics) + if seed is not None: + np.random.set_state(random_state) + toc = time.time() + self.time_evaluate = toc-tic + if "time_fit" in metrics: + test_results["time_fit"] = self.cc.time_fit + if "time_evaluate" in metrics: + test_results["time_evaluate"] = self.time_evaluate + return test_results + +class WrapRegressor(): + """ + A learner wrapped with a :class:`.ConformalRegressor` + or :class:`.ConformalPredictiveSystem`. + """ + + def __init__(self, learner): + self.cr = None + self.cps = None + self.calibrated = False + self.learner = learner + self.de = None + self.mc = None + self.seed = None + + def __repr__(self): + if self.calibrated: + if self.cr is not None: + return (f"WrapRegressor(learner={self.learner}, " + f"calibrated={self.calibrated}, " + f"predictor={self.cr})") + else: + return (f"WrapRegressor(learner={self.learner}, " + f"calibrated={self.calibrated}, " + f"predictor={self.cps})") + else: + return f"WrapRegressor(learner={self.learner}, calibrated={self.calibrated})" + + def fit(self, X, y, **kwargs): + """ + Fit learner. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features), + set of objects + y : array-like of shape (n_samples,), + target values + kwargs : optional arguments + any additional arguments are forwarded to the + ``fit`` method of the ``learner`` object + + Returns + ------- + None + + Examples + -------- + Assuming ``X_train`` and ``y_train`` to be an array and vector + with training objects and labels, respectively, a random + forest may be wrapped and fitted by: .. code-block:: python - p_values = w.predict_cps(X_test, y=37) + from sklearn.ensemble import RandomForestRegressor + from crepes import WrapRegressor - The 90th and 95th percentiles can be obtained by: + rf = WrapRegressor(RandomForestRegressor()) + rf.fit(X_train, y_train) + + Note + ---- + The learner, which can be accessed by ``rf.learner``, may be fitted + before as well as after being wrapped. + + Note + ---- + All arguments, including any additional keyword arguments, to + :meth:`.fit` are forwarded to the ``fit`` method of the learner. + """ + if isinstance(y, pd.Series): + y = y.values + self.learner.fit(X, y, **kwargs) + + def predict(self, X): + """ + Predict with learner. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features), + set of objects + + Returns + ------- + y : array-like of shape (n_samples,), + values predicted using the ``predict`` + method of the ``learner`` object. + + Examples + -------- + Assuming ``w`` is a :class:`.WrapRegressor` object for which the wrapped + learner ``w.learner`` has been fitted, (point) predictions of the + learner can be obtained for a set of test objects ``X_test`` by: .. code-block:: python - percentiles = w.predict_cps(X_test, higher_percentiles=[90,95]) + y_hat = w.predict(X_test) + + The above is equivalent to: - In the above example, the nearest higher value is returned, if there is - no value that corresponds exactly to the requested percentile. If we - instead would like to retrieve the nearest lower value, we should - write: + .. code-block:: python + + y_hat = w.learner.predict(X_test) + """ + return self.learner.predict(X) + + def calibrate(self, X, y, de=None, mc=None, oob=False, cps=False, + seed=None): + """ + Fit a :class:`.ConformalRegressor` or + :class:`.ConformalPredictiveSystem` using the wrapped learner. + + Parameters + ---------- + X : array-like of shape (n_samples, n_features), + set of objects + y : array-like of shape (n_samples,), + target values + de: :class:`crepes.extras.DifficultyEstimator`, default=None + object used for computing difficulty estimates + mc: function or :class:`crepes.extras.MondrianCategorizer`, default=None + function or :class:`crepes.extras.MondrianCategorizer` for computing Mondrian categories + oob : bool, default=False + use out-of-bag estimation + cps : bool, default=False + if cps=False, the method fits a :class:`.ConformalRegressor` + and otherwise, a :class:`.ConformalPredictiveSystem` + seed : int, default=None + set random seed + + Returns + ------- + self : object + The :class:`.WrapRegressor` object is updated with a fitted + :class:`.ConformalRegressor` or :class:`.ConformalPredictiveSystem` + + Examples + -------- + Assuming ``X_cal`` and ``y_cal`` to be an array and vector, + respectively, with objects and labels for the calibration set, + and ``w`` is a :class:`.WrapRegressor` object for which the learner + has been fitted, a standard conformal regressor is formed by: .. code-block:: python - percentiles = w.predict_cps(X_test, lower_percentiles=[90,95]) + w.calibrate(X_cal, y_cal) - The following returns prediction intervals at the 95% confidence level, - where the intervals are lower-bounded by 0: + Assuming that ``de`` is a fitted difficulty estimator, + a normalized conformal regressor is obtained by: .. code-block:: python - intervals = w.predict_cps(X_test, - lower_percentiles=2.5, - higher_percentiles=97.5, - y_min=0) + w.calibrate(X_cal, y_cal, de=de) - If we would like to obtain the conformal distributions, we could write - the following: + Assuming that ``get_categories`` is a function that returns categories + (bin labels), a Mondrian conformal regressor is obtained by: .. code-block:: python - cpds = w.predict_cps(X_test, return_cpds=True) + w.calibrate(X_cal, y_cal, mc=get_categories) - The output of the above will be an array with a row for each test - instance and a column for each calibration instance (residual). - If the learner is wrapped with a Mondrian conformal predictive system, - the above will instead result in a vector, in which each element is a - vector, as the number of calibration instances may vary between - categories. If we instead would like an array for each category, this - can be obtained by: + A normalized Mondrian conformal regressor is generated in the + following way: .. code-block:: python - cpds = w.predict_cps(X_test, return_cpds=True, cpds_by_bins=True) + w.calibrate(X_cal, y_cal, de=de, mc=get_categories) - Note - ---- - This method is available only if the learner has been wrapped with a - :class:`.ConformalPredictiveSystem`, i.e., :meth:`.calibrate` - has been called with the option ``cps=True``. + By providing the option ``oob=True``, the conformal regressor + will be calibrating using out-of-bag predictions, allowing + the full set of training objects (``X_train``) and labels (``y_train``) + to be used, e.g., + + .. code-block:: python + + w.calibrate(X_train, y_train, oob=True) + + By providing the option ``cps=True``, each of the above calls will instead + generate a :class:`.ConformalPredictiveSystem`, e.g., + + .. code-block:: python + + w.calibrate(X_cal, y_cal, de=de, cps=True) Note ---- - In case the calibration set is too small for the specified lower and - higher percentiles, a warning will be issued and the output will be - ``y_min`` and ``y_max``, respectively. - + By providing a random seed, e.g., ``seed=123``, the call to ``calibrate`` + as well as calls to the methods ``predict_int``, ``predict_cps`` and + ``evaluate`` of the :class:`.WrapRegressor` object will be deterministic. + Note ---- - Setting ``return_cpds=True`` may consume a lot of memory, as a matrix is - generated for which the number of elements is the product of the number - of calibration and test objects, unless a Mondrian approach is employed; - for the latter, this number is reduced by increasing the number of bins. + Enabling out-of-bag calibration, i.e., setting ``oob=True``, requires + that the wrapped learner has an attribute ``oob_prediction_``, which + e.g., is the case for a ``sklearn.ensemble.RandomForestRegressor``, if + enabled when created, e.g., ``RandomForestRegressor(oob_score=True)`` Note ---- - Setting ``cpds_by_bins=True`` has an effect only for Mondrian conformal - predictive systems. + The use of out-of-bag calibration, as enabled by ``oob=True``, + does not come with the theoretical validity guarantees of the regular + (inductive) conformal regressors and predictive systems, due to that + calibration and test instances are not handled in exactly the same way. """ + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + self.seed = seed if isinstance(y, pd.Series): y = y.values - if self.cps is None: - raise RuntimeError(("predict_cps requires that calibrate has been" - "called first with cps=True")) + if oob: + residuals = y - self.learner.oob_prediction_ else: - y_hat = self.learner.predict(X) - if self.de is None: - sigmas = None - else: - sigmas = self.de.apply(X) - if self.mc is None: - bins = None - elif isinstance(self.mc, MondrianCategorizer): - bins = self.mc.apply(X) - else: - bins = self.mc(X) - return self.cps.predict(y_hat, sigmas=sigmas, bins=bins, - y=y, lower_percentiles=lower_percentiles, - higher_percentiles=higher_percentiles, - y_min=y_min, y_max=y_max, - return_cpds=return_cpds, - cpds_by_bins=cpds_by_bins) - - def evaluate(self, X, y, confidence=0.95, - y_min=-np.inf, y_max=np.inf, metrics=None): + residuals = y - self.predict(X) + if de is None: + sigmas = None + else: + sigmas = de.apply(X) + self.de = de + if mc is None: + bins = None + elif isinstance(mc, MondrianCategorizer): + bins = mc.apply(X) + else: + bins = mc(X) + self.mc = mc + if not cps: + self.cr = ConformalRegressor() + self.cr.fit(residuals, sigmas=sigmas, bins=bins) + self.cps = None + else: + self.cps = ConformalPredictiveSystem() + self.cps.fit(residuals, sigmas=sigmas, bins=bins) + self.cr = None + self.calibrated = True + if seed is not None: + np.random.set_state(random_state) + return self + + def predict_int(self, X, confidence=0.95, y_min=-np.inf, y_max=np.inf, + seed=None): """ - Evaluate :class:`.ConformalRegressor` or + Predict interval using fitted :class:`.ConformalRegressor` or :class:`.ConformalPredictiveSystem`. Parameters ---------- - X : array-like of shape (n_samples, n_features) + X : array-like of shape (n_samples, n_features), set of objects - y : array-like of shape (n_samples,) - correct target values confidence : float in range (0,1), default=0.95 confidence level y_min : float or int, default=-numpy.inf minimum value to include in prediction intervals y_max : float or int, default=numpy.inf maximum value to include in prediction intervals - metrics : a string or a list of strings, default=list of all - metrics; for a learner wrapped with a conformal regressor - these are "error", "eff_mean","eff_med", "time_fit", and - "time_evaluate", while if wrapped with a conformal predictive - system, the metrics also include "CRPS". - + seed : int, default=None + set random seed + Returns ------- - results : dictionary with a key for each selected metric - estimated performance using the metrics + intervals : ndarray of shape (n_samples, 2) + prediction intervals Examples -------- - Assuming that ``X_test`` is a set of test objects, ``y_test`` is a - vector with true targets, and ``w`` is a calibrated - :class:`.WrapRegressor` object, then the latter can be evaluated at - the 90% confidence level with respect to error, mean and median - efficiency (interval size) by: + Assuming that ``X_test`` is a set of test objects and ``w`` is a + :class:`.WrapRegressor` object that has been calibrated, i.e., + :meth:`.calibrate` has been applied, prediction intervals at the + 99% confidence level can be obtained by: .. code-block:: python - results = w.evaluate(X_test, y_test, confidence=0.9, - metrics=["error", "eff_mean", "eff_med"]) + intervals = w.predict_int(X_test, confidence=0.99) - Note - ---- - If included in the list of metrics, "CRPS" (continuous-ranked - probability score) will be ignored if the :class:`.WrapRegressor` object - has been calibrated with the (default) option ``cps=False``, i.e., the - learner is wrapped with a :class:`.ConformalRegressor`. + The following provides prediction intervals at the default confidence + level (95%), where the intervals are lower-bounded by 0: + + .. code-block:: python + + intervals = w.predict_int(X_test, y_min=0) Note ---- - The use of the metric ``CRPS`` may consume a lot of memory, as a matrix - is generated for which the number of elements is the product of the - number of calibration and test objects, unless a Mondrian approach is - employed; for the latter, this number is reduced by increasing the number - of categories. + In case the specified confidence level is too high in relation to the + size of the calibration set, a warning will be issued and the output + intervals will be of maximum size. Note ---- - The reported result for ``time_fit`` only considers fitting the - conformal regressor or predictive system; not for fitting the - learner. + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. """ - if isinstance(y, pd.Series): - y = y.values if not self.calibrated: - raise RuntimeError(("evaluate requires that calibrate has been" + raise RuntimeError(("predict_int requires that calibrate has been" "called first")) else: + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) y_hat = self.learner.predict(X) if self.de is None: sigmas = None @@ -1905,363 +2268,200 @@ def evaluate(self, X, y, confidence=0.95, else: bins = self.mc(X) if self.cr is not None: - return self.cr.evaluate(y_hat, y, sigmas=sigmas, - bins=bins, confidence=confidence, - y_min=y_min, y_max=y_max) - else: - return self.cps.evaluate(y_hat, y, sigmas=sigmas, - bins=bins, confidence=confidence, + result = self.cr.predict(y_hat, sigmas=sigmas, bins=bins, + confidence=confidence, y_min=y_min, y_max=y_max) + else: + lower_percentile = (1-confidence)/2*100 + higher_percentile = (confidence+(1-confidence)/2)*100 + result = self.cps.predict(y_hat, sigmas=sigmas, bins=bins, + lower_percentiles=lower_percentile, + higher_percentiles=higher_percentile, + y_min=y_min, y_max=y_max) + if seed is not None: + np.random.set_state(random_state) + return result -class WrapClassifier(): - """ - A learner wrapped with a :class:`.ConformalClassifier`. - """ - - def __init__(self, learner): - self.cc = None - self.nc = None - self.calibrated = False - self.learner = learner - - def __repr__(self): - if self.calibrated: - return (f"WrapClassifier(learner={self.learner}, " - f"calibrated={self.calibrated}, " - f"predictor={self.cc})") - else: - return f"WrapClassifier(learner={self.learner}, calibrated={self.calibrated})" - - def fit(self, X, y, **kwargs): + def predict_cps(self, X, y=None, lower_percentiles=None, + higher_percentiles=None, y_min=-np.inf, y_max=np.inf, + return_cpds=False, cpds_by_bins=False, seed=None): """ - Fit learner. + Predict using :class:`.ConformalPredictiveSystem`. Parameters ---------- X : array-like of shape (n_samples, n_features), set of objects - y : array-like of shape (n_samples,), - target values - kwargs : optional arguments - any additional arguments are forwarded to the - ``fit`` method of the ``learner`` object + y : float, int or array-like of shape (n_samples,), default=None + values for which p-values should be returned + lower_percentiles : array-like of shape (l_values,), default=None + percentiles for which a lower value will be output + in case a percentile lies between two values + (similar to `interpolation="lower"` in `numpy.percentile`) + higher_percentiles : array-like of shape (h_values,), default=None + percentiles for which a higher value will be output + in case a percentile lies between two values + (similar to `interpolation="higher"` in `numpy.percentile`) + y_min : float or int, default=-numpy.inf + The minimum value to include in prediction intervals. + y_max : float or int, default=numpy.inf + The maximum value to include in prediction intervals. + return_cpds : Boolean, default=False + specifies whether conformal predictive distributions (cpds) + should be output or not + cpds_by_bins : Boolean, default=False + specifies whether the output cpds should be grouped by bin or not; + only applicable when bins is not None and return_cpds = True + seed : int, default=None + set random seed Returns ------- - None + results : ndarray of shape (n_samples, n_cols) or (n_samples,) + the shape is (n_samples, n_cols) if n_cols > 1 and otherwise + (n_samples,), where n_cols = p_values+l_values+h_values where + p_values = 1 if y is not None and 0 otherwise, l_values are the + number of lower percentiles, and h_values are the number of higher + percentiles. Only returned if n_cols > 0. + cpds : ndarray of (n_samples, c_values), ndarray of (n_samples,) + or list of ndarrays + conformal predictive distributions. Only returned if + return_cpds == True. For non-Mondrian conformal predictive systems, + the distributions are represented by a single array, where the + number of columns (c_values) is determined by the number of + residuals of the fitted conformal predictive system. For Mondrian + conformal predictive systems, the distributions are represented by + a vector of arrays, if cpds_by_bins = False, or a list of arrays, + with one element for each Mondrian category, if cpds_by_bins = True. Examples -------- - Assuming ``X_train`` and ``y_train`` to be an array and vector - with training objects and labels, respectively, a random - forest may be wrapped and fitted by: + Assuming that ``X_test`` is a set of test objects, ``y_test`` is a + vector with true targets, ``w`` is a :class:`.WrapRegressor` object + calibrated with the option ``cps=True``, p-values for the true targets + can be obtained by: .. code-block:: python - from sklearn.ensemble import RandomForestClassifier - from crepes import WrapClassifier - - rf = Wrap(RandomForestClassifier()) - rf.fit(X_train, y_train) - - Note - ---- - The learner, which can be accessed by ``rf.learner``, may be fitted - before as well as after being wrapped. - - Note - ---- - All arguments, including any additional keyword arguments, to - :meth:`.fit` are forwarded to the ``fit`` method of the learner. - """ - if isinstance(y, pd.Series): - y = y.values - self.learner.fit(X, y, **kwargs) - - - def predict(self, X): - """ - Predict with learner. - - Parameters - ---------- - X : array-like of shape (n_samples, n_features), - set of objects - - Returns - ------- - y : array-like of shape (n_samples,), - values predicted using the ``predict`` - method of the ``learner`` object. - - Examples - -------- - Assuming ``w`` is a :class:`.WrapClassifier` object for which the - wrapped learner ``w.learner`` has been fitted, (point) - predictions of the learner can be obtained for a set - of test objects ``X_test`` by: - - .. code-block:: python + p_values = w.predict_cps(X_test, y=y_test) - y_hat = w.predict(X_test) - - The above is equivalent to: + P-values with respect to some specific value, e.g., 37, can be + obtained by: .. code-block:: python - y_hat = w.learner.predict(X_test) - """ - return self.learner.predict(X) - - def predict_proba(self, X): - """ - Predict with learner. - - Parameters - ---------- - X : array-like of shape (n_samples, n_features), - set of objects - - Returns - ------- - y : array-like of shape (n_samples, n_classes), - predicted probabilities using the ``predict_proba`` - method of the ``learner`` object. - - Examples - -------- - Assuming ``w`` is a :class:`.WrapClassifier` object for which the - wrapped learner ``w.learner`` has been fitted, predicted - probabilities of the learner can be obtained for a set - of test objects ``X_test`` by: - - .. code-block:: python + p_values = w.predict_cps(X_test, y=37) - probabilities = w.predict_proba(X_test) - - The above is equivalent to: + The 90th and 95th percentiles can be obtained by: .. code-block:: python - probabilities = w.learner.predict_proba(X_test) - """ - return self.learner.predict_proba(X) - - def calibrate(self, X, y, oob=False, class_cond=False, nc=hinge, mc=None): - """ - Fit a :class:`.ConformalClassifier` using the wrapped learner. - - Parameters - ---------- - X : array-like of shape (n_samples, n_features), - set of objects - y : array-like of shape (n_samples,), - target values - oob : bool, default=False - use out-of-bag estimation - class_cond : bool, default=False - if class_cond=True, the method fits a Mondrian - :class:`.ConformalClassifier` using the class - labels as categories - nc : function, default = :func:`crepes.extras.hinge` - function to compute non-conformity scores - mc: function or :class:`crepes.extras.MondrianCategorizer`, default=None - function or :class:`crepes.extras.MondrianCategorizer` for computing Mondrian - categories - - Returns - ------- - self : object - Wrap object updated with a fitted :class:`.ConformalClassifier` + percentiles = w.predict_cps(X_test, higher_percentiles=[90,95]) - Examples - -------- - Assuming ``X_cal`` and ``y_cal`` to be an array and vector, - respectively, with objects and labels for the calibration set, - and ``w`` is a :class:`.WrapClassifier` object for which the learner - has been fitted, a standard conformal classifier can be formed by: + In the above example, the nearest higher value is returned, if there is + no value that corresponds exactly to the requested percentile. If we + instead would like to retrieve the nearest lower value, we should + write: .. code-block:: python - w.calibrate(X_cal, y_cal) + percentiles = w.predict_cps(X_test, lower_percentiles=[90,95]) - Assuming that ``get_categories`` is a function that returns a vector of - Mondrian categories (bin labels), a Mondrian conformal classifier can - be generated by: + The following returns prediction intervals at the 95% confidence level, + where the intervals are lower-bounded by 0: .. code-block:: python - w.calibrate(X_cal, y_cal, mc=get_categories) + intervals = w.predict_cps(X_test, + lower_percentiles=2.5, + higher_percentiles=97.5, + y_min=0) - By providing the option ``oob=True``, the conformal classifier - will be calibrating using out-of-bag predictions, allowing - the full set of training objects (``X_train``) and labels (``y_train``) - to be used, e.g., + If we would like to obtain the conformal distributions, we could write + the following: .. code-block:: python - w.calibrate(X_train, y_train, oob=True) + cpds = w.predict_cps(X_test, return_cpds=True) - By providing the option ``class_cond=True``, a Mondrian conformal classifier - will be formed using the class labels as categories, e.g., + The output of the above will be an array with a row for each test + instance and a column for each calibration instance (residual). + If the learner is wrapped with a Mondrian conformal predictive system, + the above will instead result in a vector, in which each element is a + vector, as the number of calibration instances may vary between + categories. If we instead would like an array for each category, this + can be obtained by: .. code-block:: python - w.calibrate(X_cal, y_cal, class_cond=True) + cpds = w.predict_cps(X_test, return_cpds=True, cpds_by_bins=True) Note ---- - Any Mondrian categorizer specified by the ``mc`` argument will be - ignored by :meth:`.calibrate`, if ``class_cond=True``, as the latter - implies that Mondrian categories are formed using the labels in ``y``. + This method is available only if the learner has been wrapped with a + :class:`.ConformalPredictiveSystem`, i.e., :meth:`.calibrate` + has been called with the option ``cps=True``. Note ---- - Enabling out-of-bag calibration, i.e., setting ``oob=True``, requires - that the wrapped learner has an attribute ``oob_decision_function_``, - which e.g., is the case for a ``sklearn.ensemble.RandomForestClassifier``, - if enabled when created, e.g., ``RandomForestClassifier(oob_score=True)`` + In case the calibration set is too small for the specified lower and + higher percentiles, a warning will be issued and the output will be + ``y_min`` and ``y_max``, respectively. Note ---- - The use of out-of-bag calibration, as enabled by ``oob=True``, does not - come with the theoretical validity guarantees of the regular (inductive) - conformal classifiers, due to that calibration and test instances are not - handled in exactly the same way. - """ - if isinstance(y, pd.Series): - y = y.values - self.cc = ConformalClassifier() - self.nc = nc - self.mc = mc - self.class_cond = class_cond - if oob: - alphas = nc(self.learner.oob_decision_function_, self.learner.classes_, y) - else: - alphas = nc(self.learner.predict_proba(X), self.learner.classes_, y) - if class_cond: - self.cc.fit(alphas, bins=y) - else: - if isinstance(mc, MondrianCategorizer): - bins = mc.apply(X) - self.cc.fit(alphas, bins=bins) - elif mc is not None: - bins = mc(X) - self.cc.fit(alphas, bins=bins) - else: - self.cc.fit(alphas) - self.calibrated = True - return self - - def predict_p(self, X): - """ - Obtain (smoothed) p-values using conformal classifier. - - Parameters - ---------- - X : array-like of shape (n_samples, n_features), - set of objects - - Returns - ------- - p-values : ndarray of shape (n_samples, n_classes) - p-values - - Examples - -------- - Assuming that ``X_test`` is a set of test objects and ``w`` is a - :class:`.WrapClassifier` object that has been calibrated, i.e., - :meth:`.calibrate` has been applied, the p-values for the test - objects are obtained by: - - .. code-block:: python - - p_values = w.predict_p(X_test) - """ - tic = time.time() - alphas = self.nc(self.learner.predict_proba(X)) - if self.class_cond: - p_values = np.array([ - self.cc.predict_p(alphas, - np.full(len(X), - self.learner.classes_[c]))[:, c] - for c in range(len(self.learner.classes_))]).T - else: - if isinstance(self.mc, MondrianCategorizer): - bins = self.mc.apply(X) - p_values = self.cc.predict_p(alphas, bins) - elif self.mc is not None: - bins = self.mc(X) - p_values = self.cc.predict_p(alphas, bins) - else: - p_values = self.cc.predict_p(alphas) - toc = time.time() - self.time_predict = toc-tic - return p_values - - def predict_set(self, X, confidence=0.95, smoothing=False): - """ - Obtain prediction sets using conformal classifier. - - Parameters - ---------- - X : array-like of shape (n_samples, n_features), - set of objects - confidence : float in range (0,1), default=0.95 - confidence level - smoothing : bool, default=False - use smoothed p-values - - Returns - ------- - prediction sets : ndarray of shape (n_values, n_classes) - prediction sets, where the value 1 (0) indicates - that the class label is included (excluded), i.e., - the corresponding p-value is less than 1-confidence - - Examples - -------- - Assuming that ``X_test`` is a set of test objects and ``w`` is a - :class:`.WrapClassifier` object that has been calibrated, i.e., - :meth:`.calibrate` has been applied, the prediction sets for the - test objects at the 99% confidence level are obtained by: - - .. code-block:: python + Setting ``return_cpds=True`` may consume a lot of memory, as a matrix is + generated for which the number of elements is the product of the number + of calibration and test objects, unless a Mondrian approach is employed; + for the latter, this number is reduced by increasing the number of bins. - prediction_sets = w.predict_set(X_test, confidence=0.99) + Note + ---- + Setting ``cpds_by_bins=True`` has an effect only for Mondrian conformal + predictive systems. Note ---- - Using smoothed p-values substantially increases computation time and - hardly has any effect on the predictions sets, except for when having - small calibration sets. + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. """ - tic = time.time() - alphas = self.nc(self.learner.predict_proba(X)) - if self.class_cond: - prediction_set = np.array([ - self.cc.predict_set(alphas, - np.full(len(X), - self.learner.classes_[c]), - confidence, smoothing)[:, c] - for c in range(len(self.learner.classes_))]).T + if isinstance(y, pd.Series): + y = y.values + if self.cps is None: + raise RuntimeError(("predict_cps requires that calibrate has been" + "called first with cps=True")) else: - if isinstance(self.mc, MondrianCategorizer): - bins = self.mc.apply(X) - elif self.mc is not None: - bins = self.mc(X) + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + y_hat = self.learner.predict(X) + if self.de is None: + sigmas = None else: + sigmas = self.de.apply(X) + if self.mc is None: bins = None - prediction_set = self.cc.predict_set(alphas, bins, confidence, - smoothing) - toc = time.time() - self.time_predict = toc-tic - return prediction_set + elif isinstance(self.mc, MondrianCategorizer): + bins = self.mc.apply(X) + else: + bins = self.mc(X) + result = self.cps.predict(y_hat, sigmas=sigmas, bins=bins, + y=y, lower_percentiles=lower_percentiles, + higher_percentiles=higher_percentiles, + y_min=y_min, y_max=y_max, + return_cpds=return_cpds, + cpds_by_bins=cpds_by_bins) + if seed is not None: + np.random.set_state(random_state) + return result - def evaluate(self, X, y, confidence=0.95, smoothing=False, - metrics=None): + def evaluate(self, X, y, confidence=0.95, y_min=-np.inf, y_max=np.inf, + metrics=None, seed=None): """ - Evaluate :class:`.ConformalClassifier`. + Evaluate :class:`.ConformalRegressor` or + :class:`.ConformalPredictiveSystem`. Parameters ---------- @@ -2271,36 +2471,50 @@ def evaluate(self, X, y, confidence=0.95, smoothing=False, correct target values confidence : float in range (0,1), default=0.95 confidence level - smoothing : bool, default=False - use smoothed p-values - metrics : a string or a list of strings, - default=list of all metrics, i.e., ["error", "avg_c", "one_c", - "empty", "time_fit", "time_evaluate"] + y_min : float or int, default=-numpy.inf + minimum value to include in prediction intervals + y_max : float or int, default=numpy.inf + maximum value to include in prediction intervals + metrics : a string or a list of strings, default=list of all + metrics; for a learner wrapped with a conformal regressor + these are "error", "eff_mean","eff_med", "time_fit", and + "time_evaluate", while if wrapped with a conformal predictive + system, the metrics also include "CRPS". + seed : int, default=None + set random seed Returns ------- results : dictionary with a key for each selected metric - estimated performance using the metrics, where "error" is the - fraction of prediction sets not containing the true class label, - "avg_c" is the average no. of predicted class labels, "one_c" is - the fraction of singleton prediction sets, "empty" is the fraction - of empty prediction sets, "time_fit" is the time taken to fit the - conformal classifier, and "time_evaluate" is the time taken for the - evaluation - + estimated performance using the metrics Examples -------- Assuming that ``X_test`` is a set of test objects, ``y_test`` is a vector with true targets, and ``w`` is a calibrated - :class:`.WrapClassifier` object, then the latter can be evaluated at - the 90% confidence level with respect to error, average prediction set - size and fraction of singleton predictions by: + :class:`.WrapRegressor` object, then the latter can be evaluated at + the 90% confidence level with respect to error, mean and median + efficiency (interval size) by: .. code-block:: python results = w.evaluate(X_test, y_test, confidence=0.9, - metrics=["error", "avg_c", "one_c"]) + metrics=["error", "eff_mean", "eff_med"]) + + Note + ---- + If included in the list of metrics, "CRPS" (continuous-ranked + probability score) will be ignored if the :class:`.WrapRegressor` object + has been calibrated with the (default) option ``cps=False``, i.e., the + learner is wrapped with a :class:`.ConformalRegressor`. + + Note + ---- + The use of the metric ``CRPS`` may consume a lot of memory, as a matrix + is generated for which the number of elements is the product of the + number of calibration and test objects, unless a Mondrian approach is + employed; for the latter, this number is reduced by increasing the number + of categories. Note ---- @@ -2310,9 +2524,8 @@ def evaluate(self, X, y, confidence=0.95, smoothing=False, Note ---- - Using smoothed p-values substantially increases computation time and - hardly has any effect on the results, except for when having small - calibration sets. + If a value for ``seed`` is given, it will take precedence over any ``seed`` + value given when calling ``calibrate``. """ if isinstance(y, pd.Series): y = y.values @@ -2320,17 +2533,30 @@ def evaluate(self, X, y, confidence=0.95, smoothing=False, raise RuntimeError(("evaluate requires that calibrate has been" "called first")) else: - if metrics is None: - metrics = ["error", "avg_c", "one_c", "empty", "time_fit", - "time_evaluate"] - tic = time.time() - prediction_sets = self.predict_set(X, confidence, smoothing) - test_results = get_test_results(prediction_sets, - self.learner.classes_, y, metrics) - toc = time.time() - self.time_evaluate = toc-tic - if "time_fit" in metrics: - test_results["time_fit"] = self.cc.time_fit - if "time_evaluate" in metrics: - test_results["time_evaluate"] = self.time_evaluate - return test_results + if seed is None: + seed = self.seed + if seed is not None: + random_state = np.random.get_state() + np.random.seed(seed) + y_hat = self.learner.predict(X) + if self.de is None: + sigmas = None + else: + sigmas = self.de.apply(X) + if self.mc is None: + bins = None + elif isinstance(self.mc, MondrianCategorizer): + bins = self.mc.apply(X) + else: + bins = self.mc(X) + if self.cr is not None: + result = self.cr.evaluate(y_hat, y, sigmas=sigmas, + bins=bins, confidence=confidence, + y_min=y_min, y_max=y_max) + else: + result = self.cps.evaluate(y_hat, y, sigmas=sigmas, + bins=bins, confidence=confidence, + y_min=y_min, y_max=y_max) + if seed is not None: + np.random.set_state(random_state) + return result diff --git a/src/crepes/extras.py b/src/crepes/extras.py index 300cbb8..738fb78 100644 --- a/src/crepes/extras.py +++ b/src/crepes/extras.py @@ -262,24 +262,19 @@ def fit(self, X=None, f=None, de=None, learner=None, no_bins=10, oob=False): Examples -------- - Assuming that ``X_prop_train`` is a proper training set, - then a difficulty estimator using the distances to the k - nearest neighbors can be formed in the following way - (here using the default ``k=25``): + Assuming that ``X_train`` is an array of shape (n_samples, n_features) + and ``get_values`` is a function that given ``X_train`` returns a vector + of values of shape (n_samples,), then a Mondrian categorizer can + be formed in the following way, where the boundaries for the Mondrian + categories are found by partitioning the values in the vector into five + equal-sized bins: .. code-block:: python - from crepes.extras import DifficultyEstimator - - de_knn_dist = DifficultyEstimator() - de_knn_dist.fit(X_prop_train) + from crepes.extras import MondrianCategorizer - Note - ---- - The use of out-of-bag calibration, as enabled by ``oob=True``, - does not come with the theoretical validity guarantees of the regular - (inductive) conformal regressors and predictive systems, due to that - calibration and test instances are not handled in exactly the same way. + mc = MondrianCategorizer() + mc.fit(X, f=get_values, no_bins=5) """ if f is not None: if X is not None: @@ -336,28 +331,18 @@ def apply(self, X): Examples -------- - Assuming ``de`` to be a fitted :class:`.DifficultyEstimator`, i.e., for - which :meth:`.fit` has earlier been called, difficulty estimates for a + Assuming ``mc`` to be a fitted :class:`.MondrianCategorizer`, i.e., for + which :meth:`.fit` has earlier been called, Mondrian categories for a set of objects ``X`` is obtained by: .. code-block:: python - difficulty_estimates = de.apply(X) + categories = mc.apply(X) - If ``de_oob`` is a :class:`.DifficultyEstimator` that has been fitted - with the option ``oob=True`` and a training set, then a call to - :meth:`.apply` without any objects will return the estimates for the - training set: - - .. code-block:: python - - oob_difficulty_estimates = de_oob.apply() - - For a difficulty estimator employing any of the k-nearest neighbor - approaches, the above will return an estimate for the difficulty - of each object in the training set computed using a leave-one-out - procedure, while for the variance-based approach the out-of-bag - predictions will instead be used. + Note + ---- + The array used when calling :meth:`.fit` must have the same number of + columns (``n_features``) as the array used as input to :meth:`.apply`. """ if self.f is not None: if self.bin_thresholds is None: @@ -412,8 +397,8 @@ def __repr__(self): else: return f"DifficultyEstimator(fitted={self.fitted})" - def fit(self, X=None, y=None, residuals=None, learner=None, k=25, - scaler=False, beta=0.01, oob=False): + def fit(self, X=None, f=None, y=None, residuals=None, learner=None, + k=25, scaler=False, beta=0.01, oob=False): """ Fit difficulty estimator. @@ -421,15 +406,19 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, ---------- X : array-like of shape (n_samples, n_features), default=None set of objects + f : function which given an array-like of shape (n_samples, n_features) + should return a vector of shape (n_samples,) of type int or float, + default=None + function used to compute difficulty estimates y : array-like of shape (n_samples,), default=None target values residuals : array-like of shape (n_samples,), default=None true target values - predicted values learner : an object with attribute ``learner.estimators_``, default=None an ensemble model where each model m in ``learner.estimators_`` has a - method ``m.predict`` + method ``m.predict`` (used only if f=None) k: int, default=25 - number of neighbors (used only if learner=None) + number of neighbors (used only if f=None and learner=None) scaler : bool, default=True use min-max-scaler on the difficulty estimates beta : int or float, default=0.01 @@ -501,6 +490,17 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, de_var = DifficultyEstimator() de_var.fit(X_proper_train, learner=learner_prop, scaler=True) + Difficulty estimates may also be computed by an externally defined + function. Assuming that ``diff_model`` is a fitted regression model, + for which the ``predict`` method gives estimates of the absolute + error for the objects in ``X_proper_train``, then normalized difficulty + estimates can be obtained from the following difficulty estimator: + + .. code-block:: python + + de_mod = DifficultyEstimator() + de_mod.fit(X_proper_train, f=diff_model.predict, scaler=True) + The :class:`.DifficultyEstimator` can also support the construction of conformal regressors and predictive systems that employ out-of-bag calibration. For the k-nearest neighbor approaches, the difficulty of @@ -537,6 +537,7 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, (inductive) conformal regressors and predictive systems, due to that calibration and test instances are not handled in exactly the same way. """ + self.f = f if isinstance(y, pd.Series): y = y.values self.y = y @@ -546,16 +547,9 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, self.beta = beta self.scaler = scaler self.oob = oob - if self.learner is None: - self.estimator_type = "knn" - if self.residuals is None: - if self.y is None: - self.target_type = "none" - else: - self.target_type = "labels" - else: - self.target_type = "residuals" - else: + if self.f is not None: + self.estimator_type = "function" + elif self.learner is not None: self.estimator_type = "variance" try: self.learner.estimators_ @@ -569,10 +563,53 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, raise ValueError( ("learner.estimators_ is missing the attribute " "random_state")) - if self.estimator_type == "knn": + else: + self.estimator_type = "knn" + if self.residuals is None: + if self.y is None: + self.target_type = "none" + else: + self.target_type = "labels" + else: + self.target_type = "residuals" if X is None: raise ValueError("X=None is not allowed for k-nearest" " neighbor estimators") + + if self.estimator_type == "function": + if X is None and self.scaler: + raise ValueError("X=None is allowed only if scaler=False" + " for function estimators") + if self.oob: + raise ValueError("oob=True is not allowed for function" + " estimators") + if self.scaler: + sigmas = self.f(X) + sigma_scaler = MinMaxScaler(clip=True) + sigma_scaler.fit(sigmas[:,None]) + self.sigma_scaler = sigma_scaler + + if self.estimator_type == "variance": + if X is None and (self.oob or self.scaler): + raise ValueError("X=None is allowed only if oob=False and " + "scaler=False for variance estimator") + if self.oob or self.scaler: + predictions = np.array([model.predict(X) + for model in self.learner.estimators_]) + oob_masks = np.array([ + get_oob(self.learner.estimators_[i].random_state, len(X)) + for i in range(len(self.learner.estimators_))]) + sigmas = np.array([np.var(predictions[oob_masks[:,i],i]) + for i in range(len(X))]) + if self.scaler: + sigma_scaler = MinMaxScaler(clip=True) + sigma_scaler.fit(sigmas[:,None]) + self.sigma_scaler = sigma_scaler + sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] + if self.oob: + self.sigmas = sigmas + + if self.estimator_type == "knn": nn = NearestNeighbors(n_neighbors=self.k, n_jobs=-1) nn_scaler = MinMaxScaler(clip=True) nn_scaler.fit(X) @@ -601,25 +638,7 @@ def fit(self, X=None, y=None, residuals=None, learner=None, k=25, sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] if self.oob: self.sigmas = sigmas - else: # self.estimator_type == "variance": - if X is None and (self.oob or self.scaler): - raise ValueError("X=None is allowed only if oob=False and " - "scaler=False for variance estimator") - if self.oob or self.scaler: - predictions = np.array([model.predict(X) - for model in self.learner.estimators_]) - oob_masks = np.array([ - get_oob(self.learner.estimators_[i].random_state, len(X)) - for i in range(len(self.learner.estimators_))]) - sigmas = np.array([np.var(predictions[oob_masks[:,i],i]) - for i in range(len(X))]) - if self.scaler: - sigma_scaler = MinMaxScaler(clip=True) - sigma_scaler.fit(sigmas[:,None]) - self.sigma_scaler = sigma_scaler - sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] - if self.oob: - self.sigmas = sigmas + self.fitted = True return self @@ -685,7 +704,7 @@ def apply(self, X=None): for indexes in neighbor_indexes]) if self.scaler: sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] - else: # self.estimator_type == "variance" + elif self.estimator_type == "variance": if self.oob: predictions = np.array([model.predict(X) for model in self.learner.estimators_]) @@ -700,6 +719,10 @@ def apply(self, X=None): axis=0) if self.scaler: sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] + else: # self.estimator_type == "function" + sigmas = self.f(X) + if self.scaler: + sigmas = self.sigma_scaler.transform(sigmas[:,None])[:,0] return sigmas + self.beta def get_oob(seed, n_samples):