From d04668198632a57639b972cfaddbffe167243c21 Mon Sep 17 00:00:00 2001 From: Alejandro Moreo Date: Wed, 26 Aug 2026 18:08:37 +0200 Subject: [PATCH] added bbse hard and soft and classifier adaptation wrapper --- TODO.txt | 4 +- docs/source/manuals.rst | 1 + docs/source/manuals/methods.md | 54 ++++ docs/source/quapy.classification.rst | 8 + quapy/classification/__init__.py | 1 + quapy/functional.py | 105 ++----- quapy/method/__init__.py | 6 + quapy/method/_liep_draft.py | 261 ++++++++++++++++ quapy/method/aggregative.py | 433 +++++++++++++++++++++++++-- quapy/tests/test_methods.py | 26 +- 10 files changed, 791 insertions(+), 108 deletions(-) create mode 100644 quapy/method/_liep_draft.py diff --git a/TODO.txt b/TODO.txt index 17604ab..12f1f59 100644 --- a/TODO.txt +++ b/TODO.txt @@ -19,4 +19,6 @@ scale each value by per-class thresholds, i.e., [0.33*0.1, 0.33*1, 0.33*1]/sum. - [TODO] add Friedman's method and DeBias - [TODO] check ignore warning stuff check https://docs.python.org/3/library/warnings.html#temporarily-suppressing-warnings -- [TODO] nmd and md are not selectable from qp.evaluation.evaluate as a string \ No newline at end of file +- [TODO] nmd and md are not selectable from qp.evaluation.evaluate as a string +- [TODO] add https://www.kaggle.com/datasets/Cornell-University/arxiv dataset? in the paper "Online Adaptation to + Label Distribution Shift" they say it provides a natural label shift over time \ No newline at end of file diff --git a/docs/source/manuals.rst b/docs/source/manuals.rst index 97e0628..f81ab7c 100644 --- a/docs/source/manuals.rst +++ b/docs/source/manuals.rst @@ -7,6 +7,7 @@ Manuals manuals/datasets manuals/evaluation + manuals/label-shift-adaptation manuals/methods manuals/model-selection manuals/plotting diff --git a/docs/source/manuals/methods.md b/docs/source/manuals/methods.md index b721957..5be168f 100644 --- a/docs/source/manuals/methods.md +++ b/docs/source/manuals/methods.md @@ -323,6 +323,60 @@ model.fit(*train.Xy) estim_prevalence = model.predict(test.X) ``` +RLLS computes its importance weights directly (it implements the `ImportanceWeightQuantifier` interface); see +the {ref}`Label Shift Adaptation manual ` for how to +access these weights, or use them to adapt a classifier itself rather than only estimating prevalence. + +### Black Box Shift Estimation (BBSE) + +`BBSEhard` and `BBSEsoft` are available at `qp.method.aggregative.BBSEhard` and +`qp.method.aggregative.BBSEsoft`, respectively, and implement the Black Box Shift Estimator +proposed in: + +_Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). Detecting and correcting for label shift +with black box predictors. In International conference on machine learning +(pp. 3122-3130). PMLR._ ([link to paper](https://proceedings.mlr.press/v80/lipton18a.html)) + +BBSE is similar in spirit to ACC and PACC in that it exploits the label-shift invariance of +`P(hat{Y}|Y)` to correct for the change in class prevalence between the training and the test +distributions. However, while ACC solves the linear system `q = Mp` (with `M` the matrix of +class-conditional misclassification rates and `p` the sought prevalence vector), BBSE instead +solves `q = Cw` for the importance-weight vector `w`, with `w_i = Q(i)/P(i)` the ratio between +the target and the source class priors, and `C` the joint-distribution matrix +`C_ij = P(hat{Y}=i, Y=j)` estimated on a validation split. The target prevalence estimate is +then recovered as `Q(y) = w_y * P(y)`. + +`BBSEhard` estimates `C` from crisp classifier predictions (i.e., a standard confusion matrix, +normalized to sum to 1), while `BBSEsoft` estimates it from the classifier's posterior +probabilities instead, in the same spirit in which PACC generalizes ACC. + +```python +import quapy as qp +from quapy.method.aggregative import BBSEhard, BBSEsoft +from sklearn.linear_model import LogisticRegression + +train, test = qp.datasets.fetch_UCIBinaryDataset('haberman').train_test + +model = BBSEhard(LogisticRegression(max_iter=2000), val_split=5) +model.fit(*train.Xy) +estim_prevalence = model.predict(test.X) + +# or, using posterior probabilities instead of crisp counts: +model = BBSEsoft(LogisticRegression(max_iter=2000), val_split=5) +model.fit(*train.Xy) +estim_prevalence = model.predict(test.X) +``` + +As with ACC and RLLS, both variants require validation predictions and therefore expect +`val_split` to be set whenever `fit_classifier=True`. They also accept the same `solver` +(`"minimize"`, `"exact-raise"`, `"exact-cc"`) and `norm` (`"clip"`, `"mapsimplex"`, +`"condsoftmax"`) arguments discussed above for ACC/PACC. + +Like RLLS, both `BBSEhard` and `BBSEsoft` implement the `ImportanceWeightQuantifier` interface; see the +{ref}`Label Shift Adaptation manual ` for how to +access these weights directly, or use them to adapt a classifier itself rather than only estimating +prevalence. + ### Distribution Matching Distribution Matching (DM) methods search for the mixture parameter (the sought class prevalence values) diff --git a/docs/source/quapy.classification.rst b/docs/source/quapy.classification.rst index cfc7d9b..c07b72e 100644 --- a/docs/source/quapy.classification.rst +++ b/docs/source/quapy.classification.rst @@ -12,6 +12,14 @@ quapy.classification.calibration module :undoc-members: :show-inheritance: +quapy.classification.labelshift module +--------------------------------------- + +.. automodule:: quapy.classification.labelshift + :members: + :undoc-members: + :show-inheritance: + quapy.classification.methods module ----------------------------------- diff --git a/quapy/classification/__init__.py b/quapy/classification/__init__.py index b63e0eb..33551b7 100644 --- a/quapy/classification/__init__.py +++ b/quapy/classification/__init__.py @@ -1,3 +1,4 @@ from . import calibration +from . import labelshift from . import methods from . import svmperf diff --git a/quapy/functional.py b/quapy/functional.py index 392309e..84504ee 100644 --- a/quapy/functional.py +++ b/quapy/functional.py @@ -462,7 +462,7 @@ def argmin_prevalence(loss: Callable, raise NotImplementedError() -def optim_minimize(loss: Callable, n_classes: int, return_loss=False): +def optim_minimize(loss: Callable, n_classes: int, x0='uniform', bounds='simplex', constraints='simplex', return_loss=False): """ Searches for the optimal prevalence values, i.e., an `n_classes`-dimensional vector of the (`n_classes`-1)-simplex that yields the smallest lost. This optimization is carried out by means of a constrained search using scipy's @@ -470,19 +470,32 @@ def optim_minimize(loss: Callable, n_classes: int, return_loss=False): :param loss: (callable) the function to minimize :param n_classes: (int) the number of classes, i.e., the dimensionality of the prevalence vector + :param x0: initial solution; if the string 'uniform' is passed (default) then the initial solution is the + uniform distribution; otherwise, a valid object must be provided + :param bounds: the bounds of the search space; if the string `simplex` is passed (default) then the bounds + of a simplex of appropriate dimension is instantiated; otherwise, a valid tuple must be provided + :param constraints: the constraints of valid solutions; if the string `simplex` is passed (default) then the + constraint that a point must lie on the simplex is assumed; otherwise, valid constrains must be provided + (see scipy.optimize) :param return_loss: bool, if True, returns also the value of the loss (default is False). :return: (ndarray) the best prevalence vector found or a tuple which also contains the value of the loss if return_loss=True """ from scipy import optimize - # the initial point is set as the uniform distribution - uniform_distribution = uniform_prevalence(n_classes=n_classes) + if isinstance(x0, str) and x0=='uniform': + # the initial point is set as the uniform distribution + x0 = uniform_prevalence(n_classes=n_classes) - # solutions are bounded to those contained in the unit-simplex - bounds = tuple((0, 1) for _ in range(n_classes)) # values in [0,1] - constraints = ({'type': 'eq', 'fun': lambda x: 1 - sum(x)}) # values summing up to 1 - r = optimize.minimize(loss, x0=uniform_distribution, method='SLSQP', bounds=bounds, constraints=constraints) + if isinstance(bounds, str) and bounds=='simplex': + # solutions are bounded to those contained in the unit-simplex + bounds = tuple((0, 1) for _ in range(n_classes)) # values in [0,1] + + if isinstance(constraints, str) and constraints=='simplex': + # the point must lie on the simplex + constraints = ({'type': 'eq', 'fun': lambda x: 1 - sum(x)}) # values summing up to 1 + + r = optimize.minimize(loss, x0=x0, method='SLSQP', bounds=bounds, constraints=constraints) if return_loss: return r.x, r.fun @@ -623,84 +636,6 @@ def solve_adjustment_binary(prevalence_estim: ArrayLike, tpr: float, fpr: float, return adjusted -def solve_adjustment( - class_conditional_rates: np.ndarray, - unadjusted_counts: np.ndarray, - method: Literal["inversion", "invariant-ratio"], - solver: Literal["exact", "minimize", "exact-raise", "exact-cc"]) -> np.ndarray: - """ - Function that tries to solve for :math:`p` the equation :math:`q = M p`, where :math:`q` is the vector of - `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of - :math:`P(\hat{Y}=y_i)`, and where :math:`M` is the matrix of `class-conditional rates` with :math:`M_{ij}` an - estimate of :math:`P(\hat{Y}=y_i|Y=y_j)`. - - :param class_conditional_rates: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate - of :math:`P(\hat{Y}=y_i|Y=y_j)`, that is, the probability that an instance that belongs to class :math:`y_j` - ends up being classified as belonging to class :math:`y_i` - - :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as - estimated by CC or PCC) - - :param str method: indicates the adjustment method to be used. Valid options are: - - * `inversion`: tries to solve the equation :math:`q = M p` as :math:`p = M^{-1} q` where - :math:`M^{-1}` is the matrix inversion of :math:`M`. This inversion may not exist in - degenerated cases. - * `invariant-ratio`: invariant ratio estimator of `Vaz et al. 2018 `_, - which replaces the last equation in :math:`M` with the normalization condition (i.e., that the sum of - all prevalence values must equal 1). - - :param str solver: the method to use for solving the system of linear equations. Valid options are: - - * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank - strictly lower than `n_classes`. - * `exact-cc`: if the matrix is not full rank, returns :math:`q` (i.e., the unadjusted counts) as the estimates - * `exact`: deprecated, defaults to 'exact-cc' (will be removed in future versions) - * `minimize`: minimizes a loss, so the solution always exists - """ - if solver == "exact": - warnings.warn( - "The 'exact' solver is deprecated. Use 'exact-raise' or 'exact-cc'", DeprecationWarning, stacklevel=2) - solver = "exact-cc" - - A = np.asarray(class_conditional_rates, dtype=float) - B = np.asarray(unadjusted_counts, dtype=float) - - if method == "inversion": - pass # We leave A and B unchanged - elif method == "invariant-ratio": - # Change the last equation to replace it with the normalization condition; - # copy first so this does not mutate the caller's arrays (np.asarray above - # returns the same object, not a copy, when the input is already float64) - A = A.copy() - B = B.copy() - A[-1, :] = 1.0 - B[-1] = 1.0 - else: - raise ValueError(f"unknown {method=}") - - if solver == "minimize": - def loss(prev): - return np.linalg.norm(A @ prev - B) - return optim_minimize(loss, n_classes=A.shape[0]) - elif solver in ["exact-raise", "exact-cc"]: - # Solvers based on matrix inversion, so we use try/except block - try: - return np.linalg.solve(A, B) - except np.linalg.LinAlgError: - # The matrix is not invertible. - # Depending on the solver, we either raise an error - # or return the classifier predictions without adjustment - if solver == "exact-raise": - raise - elif solver == "exact-cc": - return unadjusted_counts - else: - raise ValueError(f"Solver {solver} not known.") - else: - raise ValueError(f'unknown {solver=}') - - # ------------------------------------------------------------------------------------------ # Transformations from Compositional analysis # ------------------------------------------------------------------------------------------ diff --git a/quapy/method/__init__.py b/quapy/method/__init__.py index 0bcab12..695ac5f 100644 --- a/quapy/method/__init__.py +++ b/quapy/method/__init__.py @@ -15,6 +15,8 @@ AGGREGATIVE_METHODS = { aggregative.ACC, aggregative.PCC, aggregative.PACC, + aggregative.BBSEhard, + aggregative.BBSEsoft, aggregative.RLLS, aggregative.EMQ, aggregative.HDy, @@ -30,6 +32,7 @@ AGGREGATIVE_METHODS = { aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, + aggregative.LEIP, # aggregative.OneVsAllAggregative, confidence.BayesianCC, _bayesian.BayesianKDEy, @@ -54,12 +57,15 @@ MULTICLASS_METHODS = { aggregative.ACC, aggregative.PCC, aggregative.PACC, + aggregative.BBSEhard, + aggregative.BBSEsoft, aggregative.RLLS, aggregative.EMQ, aggregative.EDy, aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, + aggregative.LEIP, confidence.BayesianCC, _bayesian.BayesianKDEy, _bayesian.BayesianMAPLS, diff --git a/quapy/method/_liep_draft.py b/quapy/method/_liep_draft.py new file mode 100644 index 0000000..62c9bf7 --- /dev/null +++ b/quapy/method/_liep_draft.py @@ -0,0 +1,261 @@ +import numpy as np + + +def normalize(v, eps=1e-12): + v = np.asarray(v, dtype=float) + v = np.maximum(v, eps) + return v / v.sum() + + +def bayes_prior_update(probs, new_prior, source_prior, eps=1e-12): + """ + Applies the standard label-shift prior correction: + + p_new(y|x) ∝ p_old(y|x) * new_prior(y) / source_prior(y) + + Parameters + ---------- + probs : array, shape (n_samples, n_classes) or (n_classes,) + Probabilistic outputs of the classifier. + new_prior : array, shape (n_classes,) + Current estimate of the target prior. + source_prior : array, shape (n_classes,) + Source/training prior. + """ + probs = np.asarray(probs, dtype=float) + one_dim = probs.ndim == 1 + if one_dim: + probs = probs[None, :] + + new_prior = normalize(new_prior, eps=eps) + source_prior = normalize(source_prior, eps=eps) + + weights = new_prior / np.maximum(source_prior, eps) + updated = probs * weights[None, :] + updated = updated / np.maximum(updated.sum(axis=1, keepdims=True), eps) + + return updated[0] if one_dim else updated + + +def confusion_statistic(y_true, y_pred, n_classes, mode="recall", eps=1e-12): + """ + Computes the minimum diagonal statistic used to choose tau. + + mode="recall": + diag(C) / row sums, i.e. per-true-class recall. + + mode="precision": + diag(C) / column sums, i.e. per-predicted-class precision. + + The paper calls this quantity 'minimum recall', but also describes + a column-normalized confusion matrix, which would correspond more + closely to precision. I expose both options. + """ + y_true = np.asarray(y_true, dtype=int) + y_pred = np.asarray(y_pred, dtype=int) + + C = np.zeros((n_classes, n_classes), dtype=float) + for yt, yp in zip(y_true, y_pred): + C[yt, yp] += 1.0 + + if mode == "recall": + denom = C.sum(axis=1) + elif mode == "precision": + denom = C.sum(axis=0) + else: + raise ValueError("mode must be 'recall' or 'precision'.") + + diag = np.diag(C) + valid = denom > 0 + + if not np.any(valid): + raise ValueError("No valid classes found in the confusion matrix.") + + scores = diag[valid] / np.maximum(denom[valid], eps) + return float(np.min(scores)) + + +def choose_tau_from_validation( + target_probs, + val_probs, + y_val, + mode="recall", + retain_fraction=None, +): + """ + Chooses tau following the spirit of Section 4.1. + + The paper says tau is selected from the top n percentile of target + confidences, with n = min-recall * 100. To make the selected set A + contain approximately min-recall fraction of the target samples, we set: + + tau = quantile(max_probs, 1 - min_recall) + + so that about min_recall of the target points satisfy max_prob >= tau. + + If retain_fraction is provided, it overrides the validation-derived value. + """ + target_probs = np.asarray(target_probs, dtype=float) + val_probs = np.asarray(val_probs, dtype=float) + y_val = np.asarray(y_val, dtype=int) + + n_classes = target_probs.shape[1] + y_val_pred = val_probs.argmax(axis=1) + + if retain_fraction is None: + retain_fraction = confusion_statistic( + y_true=y_val, + y_pred=y_val_pred, + n_classes=n_classes, + mode=mode, + ) + + retain_fraction = float(np.clip(retain_fraction, 0.0, 1.0)) + + target_conf = target_probs.max(axis=1) + + if retain_fraction <= 0: + tau = np.inf + elif retain_fraction >= 1: + tau = -np.inf + else: + tau = np.quantile(target_conf, 1.0 - retain_fraction) + + return tau, retain_fraction + + +def leip( + target_probs, + source_prior, + tau=None, + val_probs=None, + y_val=None, + threshold_mode="recall", + count_smoothing=0.0, + eps=1e-12, + return_details=False, +): + """ + LEIP: Label shift Estimation with Incremental Prior update. + + Parameters + ---------- + target_probs : array, shape (n_target, n_classes) + Probabilistic classifier outputs on the target/test set. + source_prior : array, shape (n_classes,) + Source class prior p_s(y). + tau : float or None + Confidence threshold. If None, it is estimated using validation data. + val_probs : array, shape (n_val, n_classes), optional + Validation probabilistic outputs, required if tau is None. + y_val : array, shape (n_val,), optional + Validation labels, required if tau is None. + threshold_mode : {"recall", "precision"} + Statistic used to choose tau from validation data. + count_smoothing : float + Optional additive smoothing for pseudo-label counts. Set to 0.0 + for a closer implementation of the paper; use a small value such + as 1e-8 for extra numerical robustness. + eps : float + Numerical stabilizer. + return_details : bool + If True, returns diagnostic information. + + Returns + ------- + estimated_prior : array, shape (n_classes,) + Estimated target class distribution. + details : dict, optional + Returned only if return_details=True. + """ + target_probs = np.asarray(target_probs, dtype=float) + if target_probs.ndim != 2: + raise ValueError("target_probs must have shape (n_samples, n_classes).") + + n_target, n_classes = target_probs.shape + source_prior = normalize(source_prior, eps=eps) + + if n_classes != len(source_prior): + raise ValueError("source_prior must have one entry per class.") + + # Step 1: choose tau if needed + if tau is None: + if val_probs is None or y_val is None: + raise ValueError("val_probs and y_val are required when tau is None.") + + tau, retain_fraction = choose_tau_from_validation( + target_probs=target_probs, + val_probs=val_probs, + y_val=y_val, + mode=threshold_mode, + ) + else: + retain_fraction = None + + target_conf = target_probs.max(axis=1) + target_top = target_probs.argmax(axis=1) + + # Step 2: high-confidence set A + A_mask = target_conf >= tau + A_labels = target_top[A_mask] + + counts = np.full(n_classes, count_smoothing, dtype=float) + + if len(A_labels) > 0: + counts += np.bincount(A_labels, minlength=n_classes) + current_prior = counts / counts.sum() + else: + # Fallback if tau is too strict. + # One could also use classify-and-count over the full target set. + current_prior = source_prior.copy() + + # Step 3: low-confidence set B, sorted by decreasing confidence + B_indices = np.where(~A_mask)[0] + B_indices = B_indices[np.argsort(-target_conf[B_indices])] + + # Step 4: incremental pass over B + incremental_labels = [] + + for idx in B_indices: + corrected = bayes_prior_update( + probs=target_probs[idx], + new_prior=current_prior, + source_prior=source_prior, + eps=eps, + ) + pseudo_label = int(np.argmax(corrected)) + incremental_labels.append(pseudo_label) + + counts[pseudo_label] += 1.0 + current_prior = counts / counts.sum() + + estimated_intermediate_prior = current_prior.copy() + + # Step 5: final complete pass over all target instances + corrected_all = bayes_prior_update( + probs=target_probs, + new_prior=estimated_intermediate_prior, + source_prior=source_prior, + eps=eps, + ) + + final_labels = corrected_all.argmax(axis=1) + estimated_prior = np.bincount(final_labels, minlength=n_classes).astype(float) + estimated_prior /= estimated_prior.sum() + + if not return_details: + return estimated_prior + + details = { + "tau": tau, + "retain_fraction": retain_fraction, + "n_A": int(A_mask.sum()), + "n_B": int((~A_mask).sum()), + "A_mask": A_mask, + "intermediate_prior": estimated_intermediate_prior, + "final_labels": final_labels, + "corrected_probs": corrected_all, + "incremental_labels": np.asarray(incremental_labels, dtype=int), + } + + return estimated_prior, details \ No newline at end of file diff --git a/quapy/method/aggregative.py b/quapy/method/aggregative.py index f678546..23341b6 100644 --- a/quapy/method/aggregative.py +++ b/quapy/method/aggregative.py @@ -368,6 +368,65 @@ class BinaryAggregativeQuantifier(AggregativeQuantifier, BinaryQuantifier): return super().fit(X, y) +class ImportanceWeightQuantifier(AggregativeQuantifier, ABC): + """ + Abstract mixin for aggregative quantifiers that estimate the target prevalence by first computing a vector + of importance weights :math:`w_y=Q(y)/P(y)` (with :math:`P` and :math:`Q` the training and target + distributions), and then rescaling the training prevalence by these weights, i.e., + :math:`\\hat{p}(y) \\propto w_y \\cdot P(y)`. + + Subclasses must set a fitted attribute `self.train_prevalence_` (typically in :meth:`aggregation_fit`) and + must implement :meth:`_weights_from_predictions`. This class provides a template implementation of + :meth:`aggregate`, together with :meth:`get_importance_weights` and :meth:`quantify_and_weigh`, none of + which mutate any internal state, so that they are all safe to call concurrently (e.g., from different + threads) on the same fitted instance for different batches of target instances. + """ + + @abstractmethod + def _weights_from_predictions(self, classif_predictions) -> np.ndarray: + """ + Computes the vector of importance weights from the classifier predictions on a batch of (unlabelled) + target instances. + + :param classif_predictions: array-like with the classifier predictions (crisp or soft, depending on + the subclass) for the target instances + :return: np.ndarray of shape `(n_classes,)` + """ + ... + + def _prevalence_from_weights(self, weights: np.ndarray) -> np.ndarray: + estimate = self.train_prevalence_ * weights + return F.normalize_prevalence(estimate, method=self.norm) + + def aggregate(self, classif_predictions): + weights = self._weights_from_predictions(classif_predictions) + return self._prevalence_from_weights(weights) + + def get_importance_weights(self, instances) -> np.ndarray: + """ + Estimates the vector of importance weights :math:`w_y=Q(y)/P(y)` for the given (unlabelled) target + instances. + + :param instances: array-like of shape `(n_instances, n_dimensions)`, the target instances + :return: np.ndarray of shape `(n_classes,)` + """ + classif_predictions = self.classify(instances) + return self._weights_from_predictions(classif_predictions) + + def quantify_and_weigh(self, instances): + """ + Jointly returns the estimated target prevalence and the importance weights used to obtain it, computed + from a single pass of classifier predictions over the given (unlabelled) target instances. + + :param instances: array-like of shape `(n_instances, n_dimensions)`, the target instances + :return: a tuple `(prevalence, weights)`, both np.ndarray of shape `(n_classes,)` + """ + classif_predictions = self.classify(instances) + weights = self._weights_from_predictions(classif_predictions) + prevalence = self._prevalence_from_weights(weights) + return prevalence, weights + + # Methods # ------------------------------------ class CC(AggregativeCrispQuantifier): @@ -456,8 +515,10 @@ class ACC(AggregativeCrispQuantifier): :param str method: adjustment method to be used: - * 'inversion': matrix inversion method based on the matrix equality :math:`P(C)=P(C|Y)P(Y)`, - which tries to invert :math:`P(C|Y)` matrix. + * 'inversion': matrix inversion method. Based on the matrix equality :math:`q=M p`, with + :math:`q` the prevalence vector estimated by CC, :math:`M` the matrix with entries :math:`i,j` representing + :math:`P(\\hat{Y}=i|Y=j)`, and :math:`p` the sought class prevalence vector, the matrix inversion + tries to solve for :math:`p=M^{-1} q` * 'invariant-ratio': invariant ratio estimator of `Vaz et al. 2018 `_, which replaces the last equation with the normalization condition. @@ -467,8 +528,7 @@ class ACC(AggregativeCrispQuantifier): strictly less than `n_classes`. * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) - * 'exact': deprecated, defaults to 'exact-cc' - * 'minimize': minimizes the L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the default parameter. More details about this can be consulted in `Bunse, M. "On Multi-Class Extensions of Adjusted Classify and Count", on proceedings of the 2nd International Workshop on Learning to Quantify: Methods and Applications (LQ 2022), ECML/PKDD 2022, Grenoble (France) @@ -489,7 +549,7 @@ class ACC(AggregativeCrispQuantifier): classifier: BaseEstimator = None, fit_classifier = True, val_split = 5, - solver: Literal['minimize', 'exact', 'exact-raise', 'exact-cc'] = 'minimize', + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', method: Literal['inversion', 'invariant-ratio'] = 'inversion', norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', n_jobs=None, @@ -500,7 +560,7 @@ class ACC(AggregativeCrispQuantifier): self.method = method self.norm = norm - SOLVERS = ['exact', 'minimize', 'exact-raise', 'exact-cc'] + SOLVERS = ['minimize', 'exact-raise', 'exact-cc'] METHODS = ['inversion', 'invariant-ratio'] NORMALIZATIONS = ['clip', 'mapsimplex', 'condsoftmax', None] @@ -551,8 +611,8 @@ class ACC(AggregativeCrispQuantifier): @classmethod def getPteCondEstim(cls, classes, y, y_): """ - Estimate the matrix with entry (i,j) being the estimate of P(hat_yi|yj), that is, the probability that a - document that belongs to yj ends up being classified as belonging to yi + Estimate the matrix with entry (i,j) being the estimate of P(hat_yi|yj), that is, the probability that an + instance that belongs to class j ends up being classified as belonging to class i :param classes: array-like with the class names :param y: array-like with the true labels @@ -570,15 +630,234 @@ class ACC(AggregativeCrispQuantifier): return conf def aggregate(self, classif_predictions): - prevs_estim = self.cc.aggregate(classif_predictions) - estimate = F.solve_adjustment( + prevs_estim_cc = self.cc.aggregate(classif_predictions) + estimate = ACC.solve_adjustment( class_conditional_rates=self.Pte_cond_estim_, - unadjusted_counts=prevs_estim, + unadjusted_counts=prevs_estim_cc, solver=self.solver, method=self.method, ) return F.normalize_prevalence(estimate, method=self.norm) + @classmethod + def solve_adjustment(cls, + class_conditional_rates: np.ndarray, + unadjusted_counts: np.ndarray, + method: Literal["inversion", "invariant-ratio"], + solver: Literal["minimize", "exact-raise", "exact-cc"]) -> np.ndarray: + """ + Function that tries to solve for :math:`p` the equation :math:`q = M p`, where :math:`q` is the vector of + `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of + :math:`P(\\hat{Y}=y_i)`, and where :math:`M` is the matrix of `class-conditional rates` with :math:`M_{ij}` an + estimate of :math:`P(\\hat{Y}=y_i|Y=y_j)`. + + :param class_conditional_rates: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate + of :math:`P(\\hat{Y}=y_i|Y=y_j)`, that is, the probability that an instance ends up being classified as + belonging to class :math:`y_i` given it actually belonged to class :math:`y_j` + + :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as + estimated by CC or PCC) + + :param str method: indicates the adjustment method to be used. Valid options are: + + * `inversion`: tries to solve the equation :math:`q = M p` as :math:`p = M^{-1} q` where + :math:`M^{-1}` is the matrix inversion of :math:`M`. This inversion may not exist in + degenerated cases. + * `invariant-ratio`: invariant ratio estimator of `Vaz et al. 2018 `_, + which replaces the last equation in :math:`M` with the normalization condition (i.e., that the sum of + all prevalence values must equal 1). + + :param str solver: the method to use for solving the system of linear equations. Valid options are: + + * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly lower than `n_classes`. + * `exact-cc`: if the matrix is not full rank, returns :math:`q` (i.e., the unadjusted counts) as the estimates + * `minimize`: minimizes a loss, so the solution always exists + """ + + A = np.asarray(class_conditional_rates, dtype=float).copy() + B = np.asarray(unadjusted_counts, dtype=float).copy() + + if method == "inversion": + pass # leave A and B unchanged + elif method == "invariant-ratio": + # Change the last equation to replace it with the normalization condition; + # copy first so this does not mutate the caller's arrays (np.asarray above + # returns the same object, not a copy, when the input is already float64) + A[-1, :] = 1.0 + B[-1] = 1.0 + else: + raise ValueError(f"unknown {method=}") + + if solver in ["exact-raise", "exact-cc"]: + try: + return np.linalg.solve(A, B) + except np.linalg.LinAlgError: + if solver=='exact-cc': + return unadjusted_counts + else: + raise + elif solver == "minimize": + def loss(prev): + return np.linalg.norm(A @ prev - B) + return F.optim_minimize(loss, n_classes=A.shape[0], return_loss=False) + else: + raise ValueError(f"unknown {solver=}") + + +class BBSEhard(ImportanceWeightQuantifier, AggregativeCrispQuantifier): + """ + `Black Box Shift Estimator` (BBSE) hard aims at finding the importance weights :math:`w_i=Q(i)/P(i)`, + with :math:`P` and :math:`Q` the training and test distributions. + BBSE is similar in spirit to ACC, but it solves the problem :math:`q=C w`, with :math:`q` the prevalence vector + estimated by CC, :math:`C` the matrix with entries :math:`i,j` representing :math:`P(\\hat{Y}=i,Y=j)`, and + :math:`w` the sought vector of importance weights. The `hard` variant estimates these quantities using + crisp counts. + + BBSE was proposed in + `Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). + Detecting and correcting for label shift with black box predictors. + In International conference on machine learning (pp. 3122-3130). PMLR. + `_. + + :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be + the one indicated in `qp.environ['DEFAULT_CLS']` + + :param fit_classifier: whether to train the learner (default is True). Set to False if the + learner has been trained outside the quantifier. + + :param val_split: specifies the data used for generating classifier predictions. This specification + can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to + be extracted from the training set; or as an integer (default 5), indicating that the predictions + are to be generated in a `k`-fold cross-validation manner (with this integer indicating the value + for `k`); or as a tuple (X,y) defining the specific set of data to use for validation. + + :param str solver: indicates the method to use for solving the system of linear equations. Valid options are: + + * 'exact-raise': tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly less than `n_classes`. + * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to + no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) + * 'exact': deprecated, defaults to 'exact-cc' + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + default parameter. More details about this can be consulted in + `Tachet des Combes, R., Zhao, H., Wang, Y. X., & Gordon, G. J. (2020). + Domain adaptation with conditional distribution matching and generalized label shift. + Advances in Neural Information Processing Systems, 33, 19276-19289. + `_. + + :param str norm: the method to use for normalization. + + * `clip`, the values are clipped to the range [0,1] and then L1-normalized. + * `mapsimplex` projects vectors onto the probability simplex. This implementation relies on + `Mathieu Blondel's projection_simplex_sort `_ + * `condsoftmax`, applies a softmax normalization only to prevalence vectors that lie outside the simplex + + :param n_jobs: number of parallel workers + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier = True, + val_split = 5, + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', + norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.n_jobs = qp._get_njobs(n_jobs) + self.solver = solver + self.norm = norm + + def _check_init_parameters(self): + if self.solver not in ACC.SOLVERS: + raise ValueError(f"unknown solver; valid ones are {ACC.SOLVERS}") + if self.norm not in ACC.NORMALIZATIONS: + raise ValueError(f"unknown normalization; valid ones are {ACC.NORMALIZATIONS}") + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimates the misclassification rates. + :param classif_predictions: array-like with the predicted labels + :param labels: array-like with the true labels associated to each predicted label + """ + true_labels = labels + pred_labels = classif_predictions + self.cc = CC(self.classifier, fit_classifier=False) + self.confusion = BBSEhard.getConfusionJointProb(self.classifier.classes_, true_labels, pred_labels) + self.train_prevalence_ = F.prevalence_from_labels(labels, classes=self.classifier.classes_) + + @classmethod + def getConfusionJointProb(cls, classes, y, y_): + """ + Estimate the matrix with entry (i,j) being the estimate of :math:`P(\\hat{Y}=i,Y=j)` + + :param classes: array-like with the class names + :param y: array-like with the true labels + :param y_: array-like with the estimated labels + :return: np.ndarray + """ + conf = confusion_matrix(y, y_, labels=classes).T + joint_probs = conf / conf.sum() + return joint_probs + + def _weights_from_predictions(self, classif_predictions): + prevs_estim_cc = self.cc.aggregate(classif_predictions) + weights = BBSEhard.solve_importance_weights( + joint_probs=self.confusion, + unadjusted_counts=prevs_estim_cc, + solver=self.solver, + ) + return np.clip(weights, 0.0, None) + + @classmethod + def solve_importance_weights(cls, + joint_probs: np.ndarray, + unadjusted_counts: np.ndarray, + solver: Literal["minimize", "exact-raise", "exact-cc"]) -> np.ndarray: + """ + Function that tries to solve for :math:`p` the equation :math:`q = C w`, where :math:`q` is the vector of + `unadjusted counts` (as estimated, e.g., via classify and count) with :math:`q_i` an estimate of + :math:`P(\\hat{Y}=y_i)`, and where :math:`C` is the confusion matrix distribution with :math:`C_{ij}` an + estimate of :math:`P(\\hat{Y}=y_i,Y=y_j)`. + + :param joint_probs: array of shape `(n_classes, n_classes,)` with entry `(i,j)` being the estimate + of :math:`P(\\hat{Y}=y_i,Y=y_j)` + + :param unadjusted_counts: array of shape `(n_classes,)` containing the unadjusted prevalence values (e.g., as + estimated by CC or PCC) + + :param str solver: the method to use for solving the system of linear equations. Valid options are: + + * `exact-raise`: tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly lower than `n_classes`. + * `exact-cc`: if the matrix is not full rank, returns a vector of ones as the weights + * `minimize`: minimizes a loss, so the solution always exists + """ + A = np.asarray(joint_probs, dtype=float) + B = np.asarray(unadjusted_counts, dtype=float) + all_ones = np.full_like(B, fill_value=1., dtype=float) + + if solver in ["exact-raise", "exact-cc"]: + try: + return np.linalg.solve(A, B) + except np.linalg.LinAlgError: + if solver=='exact-cc': + return all_ones + else: + raise + elif solver == "minimize": + def loss(prev): + return np.linalg.norm(A @ prev - B) + + n_dims = len(all_ones) + bounds = [(0, np.inf)] * n_dims + return F.optim_minimize(loss, n_classes=A.shape[0], x0=all_ones, bounds=bounds, constraints=()) + else: + raise ValueError(f"unknown {solver=}") + + class PACC(AggregativeSoftQuantifier): """ @@ -615,7 +894,6 @@ class PACC(AggregativeSoftQuantifier): Raises an error if the matrix has rank strictly less than `n_classes`. * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) - * 'exact': deprecated, defaults to 'exact-cc' * 'minimize': minimizes the L2 norm of :math:`|Ax-B|`. This one generally works better, and is the default parameter. More details about this can be consulted in `Bunse, M. "On Multi-Class Extensions of Adjusted Classify and Count", on proceedings of the 2nd International Workshop on Learning to @@ -637,7 +915,7 @@ class PACC(AggregativeSoftQuantifier): classifier: BaseEstimator = None, fit_classifier=True, val_split=5, - solver: Literal['minimize', 'exact', 'exact-raise', 'exact-cc'] = 'minimize', + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', method: Literal['inversion', 'invariant-ratio'] = 'inversion', norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', n_jobs=None @@ -671,7 +949,7 @@ class PACC(AggregativeSoftQuantifier): def aggregate(self, classif_posteriors): prevs_estim = self.pcc.aggregate(classif_posteriors) - estimate = F.solve_adjustment( + estimate = ACC.solve_adjustment( class_conditional_rates=self.Pte_cond_estim_, unadjusted_counts=prevs_estim, solver=self.solver, @@ -693,7 +971,120 @@ class PACC(AggregativeSoftQuantifier): return confusion.T -class RLLS(AggregativeSoftQuantifier): +class BBSEsoft(ImportanceWeightQuantifier, AggregativeSoftQuantifier): + """ + `Black Box Shift Estimator` (BBSE) soft, the probabilistic variant of :class:`BBSEhard` that relies on the + posterior probabilities returned by a probabilistic classifier, instead of on crisp counts, to estimate the + joint distribution :math:`P(\\hat{Y}=i,Y=j)`. As in :class:`BBSEhard`, the sought importance weights + :math:`w_i=Q(i)/P(i)` (with :math:`P` and :math:`Q` the training and test distributions) are obtained by + solving :math:`q=C w`, with :math:`q` the (now probabilistic) prevalence vector estimated by PCC and + :math:`C` the joint-probability matrix. + + BBSE was proposed in + `Lipton, Z., Wang, Y. X., & Smola, A. (2018, July). + Detecting and correcting for label shift with black box predictors. + In International conference on machine learning (pp. 3122-3130). PMLR. + `_. + + :param classifier: a scikit-learn's BaseEstimator, or None, in which case the classifier is taken to be + the one indicated in `qp.environ['DEFAULT_CLS']` + + :param fit_classifier: whether to train the learner (default is True). Set to False if the + learner has been trained outside the quantifier. + + :param val_split: specifies the data used for generating classifier predictions. This specification + can be made as float in (0, 1) indicating the proportion of stratified held-out validation set to + be extracted from the training set; or as an integer (default 5), indicating that the predictions + are to be generated in a `k`-fold cross-validation manner (with this integer indicating the value + for `k`); or as a tuple (X,y) defining the specific set of data to use for validation. + + :param str solver: indicates the method to use for solving the system of linear equations. Valid options are: + + * 'exact-raise': tries to solve the system using matrix inversion. Raises an error if the matrix has rank + strictly less than `n_classes`. + * 'exact-cc': if the matrix is not of full rank, returns `p_c` as the estimates, which corresponds to + no adjustment (i.e., the classify and count method. See :class:`quapy.method.aggregative.CC`) + * 'minimize': minimizes the squared L2 norm of :math:`|Ax-B|`. This one generally works better, and is the + default parameter. More details about this can be consulted in + `Tachet des Combes, R., Zhao, H., Wang, Y. X., & Gordon, G. J. (2020). + Domain adaptation with conditional distribution matching and generalized label shift. + Advances in Neural Information Processing Systems, 33, 19276-19289. + `_. + + :param str norm: the method to use for normalization. + + * `clip`, the values are clipped to the range [0,1] and then L1-normalized. + * `mapsimplex` projects vectors onto the probability simplex. This implementation relies on + `Mathieu Blondel's projection_simplex_sort `_ + * `condsoftmax`, applies a softmax normalization only to prevalence vectors that lie outside the simplex + + :param n_jobs: number of parallel workers + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier = True, + val_split = 5, + solver: Literal['minimize', 'exact-raise', 'exact-cc'] = 'minimize', + norm: Literal['clip', 'mapsimplex', 'condsoftmax'] = 'clip', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.n_jobs = qp._get_njobs(n_jobs) + self.solver = solver + self.norm = norm + + def _check_init_parameters(self): + if self.solver not in ACC.SOLVERS: + raise ValueError(f"unknown solver; valid ones are {ACC.SOLVERS}") + if self.norm not in ACC.NORMALIZATIONS: + raise ValueError(f"unknown normalization; valid ones are {ACC.NORMALIZATIONS}") + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimates the joint distribution P(hat_Y,Y), using posterior probabilities in place of crisp counts. + + :param classif_predictions: array-like with posterior probabilities + :param labels: array-like with the true labels associated to each vector of posterior probabilities + """ + posteriors = classif_predictions + true_labels = labels + self.pcc = PCC(self.classifier, fit_classifier=False) + self.confusion = BBSEsoft.getConfusionJointProb(self.classifier.classes_, true_labels, posteriors) + self.train_prevalence_ = F.prevalence_from_labels(labels, classes=self.classifier.classes_) + + @classmethod + def getConfusionJointProb(cls, classes, y, posteriors): + """ + Estimate the matrix with entry (i,j) being the estimate of :math:`P(\\hat{Y}=i,Y=j)`, using the + posterior probabilities of the instances belonging to class :math:`j` in place of their hard predictions. + + :param classes: array-like with the class names + :param y: array-like with the true labels + :param posteriors: array-like of shape `(n_instances, n_classes,)` with posterior probabilities + :return: np.ndarray + """ + n_classes = len(classes) + joint_probs = np.zeros((n_classes, n_classes), dtype=float) + for j, class_ in enumerate(classes): + idx = y == class_ + if idx.any(): + joint_probs[:, j] = posteriors[idx].sum(axis=0) + joint_probs /= joint_probs.sum() + return joint_probs + + def _weights_from_predictions(self, classif_posteriors): + prevs_estim_pcc = self.pcc.aggregate(classif_posteriors) + weights = BBSEhard.solve_importance_weights( + joint_probs=self.confusion, + unadjusted_counts=prevs_estim_pcc, + solver=self.solver, + ) + return np.clip(weights, 0.0, None) + + +class RLLS(ImportanceWeightQuantifier, AggregativeSoftQuantifier): """ `Regularized Learning for Domain Adaptation under Label Shifts `_, used here as an aggregative @@ -750,7 +1141,6 @@ class RLLS(AggregativeSoftQuantifier): self.delta = delta self.clip_weights = clip_weights self.norm = norm - self.last_w_ = None def _check_init_parameters(self): _get_cvxpy() @@ -781,18 +1171,15 @@ class RLLS(AggregativeSoftQuantifier): self.pz_ = _rlls_predicted_marginal(classif_predictions, mode=self.mode) self.rho_ = _rlls_compute_3deltaC(len(self.classes_), len(labels), self.delta) - def aggregate(self, classif_posteriors): + def _weights_from_predictions(self, classif_posteriors): qz = _rlls_predicted_marginal(classif_posteriors, mode=self.mode) - w = _rlls_compute_weights( + return _rlls_compute_weights( self.C_zy_, qz, self.pz_, rho=self.alpha * self.rho_, clip=self.clip_weights, ) - self.last_w_ = w - estimate = self.train_prevalence_ * w - return F.normalize_prevalence(estimate, method=self.norm) class EMQ(AggregativeSoftQuantifier): @@ -1831,6 +2218,10 @@ KDEyML = _kdey.KDEyML KDEyHD = _kdey.KDEyHD KDEyCS = _kdey.KDEyCS +from . import _liep + +LEIP = _liep.LEIP + # --------------------------------------------------------------- # aliases diff --git a/quapy/tests/test_methods.py b/quapy/tests/test_methods.py index dfd3b12..7a225dd 100644 --- a/quapy/tests/test_methods.py +++ b/quapy/tests/test_methods.py @@ -8,7 +8,7 @@ from sklearn.linear_model import LogisticRegression from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_METHODS from quapy.method.non_aggregative import DMx, EDx, HDx -from quapy.method.aggregative import ACC, DMy, EDy, KDEyCS, RLLS +from quapy.method.aggregative import ACC, BBSEhard, BBSEsoft, DMy, EDy, KDEyCS, LEIP, RLLS from quapy.method.meta import Ensemble from quapy.functional import check_prevalence_vector from quapy.tests._synthetic import make_dataset @@ -21,6 +21,7 @@ OPTIONAL_AGGREGATIVE_METHODS = { 'PQ', 'RLLS', 'EDy', + 'LEIP', } OPTIONAL_NON_AGGREGATIVE_METHODS = { @@ -178,6 +179,29 @@ class TestMethods(unittest.TestCase): self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_leip(self): + dataset = TestMethods.tiny_dataset_multiclass + q = LEIP(LogisticRegression(max_iter=2000), val_split=3) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + def test_leip_fixed_tau(self): + dataset = TestMethods.tiny_dataset_binary + q = LEIP(LogisticRegression(max_iter=2000), val_split=None, tau=0.6) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + + def test_bbse(self): + dataset = TestMethods.tiny_dataset_multiclass + for cls in (BBSEhard, BBSEsoft): + for solver in ('minimize', 'exact-raise', 'exact-cc'): + q = cls(LogisticRegression(max_iter=2000), val_split=3, solver=solver) + q.fit(*dataset.training.Xy) + estim_prevalences = q.predict(dataset.test.X) + self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_edy(self): try: import quadprog # noqa: F401