moreo
/
QuaPy
1
0
Fork 1
Code Issues Pull Requests Releases Wiki Activity

added bbse hard and soft and classifier adaptation wrapper

This commit is contained in:
Alejandro Moreo Fernandez 2026-08-26 18:08:37 +02:00
parent eeaa6776b1
commit d046681986
10 changed files with 791 additions and 108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
TODO.txt
View File

@ -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
- [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

1
docs/source/manuals.rst
View File

@ -7,6 +7,7 @@ Manuals
manuals/datasets
manuals/evaluation
manuals/label-shift-adaptation
manuals/methods
manuals/model-selection
manuals/plotting

54
docs/source/manuals/methods.md
View File

@ -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 <manuals/label-shift-adaptation:Label Shift Adaptation>` 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 <manuals/label-shift-adaptation:Label Shift Adaptation>` 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)

8
docs/source/quapy.classification.rst
View File

@ -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
-----------------------------------

1
quapy/classification/__init__.py
View File

@ -1,3 +1,4 @@
from . import calibration
from . import labelshift
from . import methods
from . import svmperf

105
quapy/functional.py
View File

@ -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 <https://jmlr.org/papers/v20/18-456.html>`_,
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
# ------------------------------------------------------------------------------------------

6
quapy/method/__init__.py
View File

@ -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,

261
quapy/method/_liep_draft.py Normal file
View File

@ -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

433
quapy/method/aggregative.py
View File

@ -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 <https://jmlr.org/papers/v20/18-456.html>`_,
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 <https://jmlr.org/papers/v20/18-456.html>`_,
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.
<https://proceedings.mlr.press/v80/lipton18a.html>`_.
: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.
<https://proceedings.neurips.cc/paper_files/paper/2020/hash/dfbfa7ddcfffeb581f50edcf9a0204bb-Abstract.html>`_.
: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 <https://gist.github.com/mblondel/6f3b7aaad90606b98f71>`_
* `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.
<https://proceedings.mlr.press/v80/lipton18a.html>`_.
: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.
<https://proceedings.neurips.cc/paper_files/paper/2020/hash/dfbfa7ddcfffeb581f50edcf9a0204bb-Abstract.html>`_.
: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 <https://gist.github.com/mblondel/6f3b7aaad90606b98f71>`_
* `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
<https://arxiv.org/abs/1903.09734>`_, 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

26
quapy/tests/test_methods.py
View File

@ -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