diff --git a/README.md b/README.md index fcb820d..4110376 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # QuaPy +## version 0.2.1 + QuaPy is an open source framework for quantification (a.k.a. supervised prevalence estimation, or learning to quantify) written in Python. @@ -74,6 +76,7 @@ See the [documentation](https://hlt-isti.github.io/QuaPy/manuals.html) for detai * Implementation of many popular quantification methods (Classify-&-Count and its variants, Expectation Maximization, quantification methods based on structured output learning, HDy, QuaNet, quantification ensembles, among others). +* Support for uncertainty quantification via bootstrap-based and Bayesian methods, including confidence intervals and simplex-aware confidence regions. * Versatile functionality for performing evaluation based on sampling generation protocols (e.g., APP, NPP, etc.). * Implementation of most commonly used evaluation metrics (e.g., AE, RAE, NAE, NRAE, SE, KLD, NKLD, etc.). * Datasets frequently used in quantification (textual and numeric), including: diff --git a/TODO.txt b/TODO.txt index de40ed9..f522a47 100644 --- a/TODO.txt +++ b/TODO.txt @@ -24,33 +24,9 @@ Solve the pre-trained classifier issues. An example is the coptic-codes script I work for having access to classes_; think also the case in which the precomputed outputs are already generated as in the unifying problems code. -Para quitar el labelledcollection de los métodos: - -- El follón viene por la semántica confusa de fit en agregativos, que recibe 3 parámetros: - - data: LabelledCollection, que puede ser: - - el training set si hay que entrenar el clasificador - - None si no hay que entregar el clasificador - - el validation, que entra en conflicto con val_split, si no hay que entrenar clasificador - - fit_classifier: dice si hay que entrenar el clasificador o no, y estos cambia la semántica de los otros - - val_split: que puede ser: - - un número: el número de kfcv, lo cual implica fit_classifier=True y data=todo el training set - - una fración en [0,1]: que indica la parte que usamos para validation; implica fit_classifier=True y data=train+val - - un labelled collection: el conjunto de validación específico; no implica fit_classifier=True ni False -- La forma de quitar la dependencia de los métodos con LabelledCollection debería ser así: - - En el constructor se dice si el clasificador que se recibe por parámetro hay que entrenarlo o ya está entrenado; - es decir, hay un fit_classifier=True o False. - - fit_classifier=True: - - data en fit es todo el training incluyendo el validation y todo - - val_split: - - int: número de folds en kfcv - - proporción en [0,1] - - fit_classifier=False: - - - [TODO] document confidence in manuals - [TODO] Test the return_type="index" in protocols and finish the "distributing_samples.py" example -- [TODO] Add EDy (an implementation is available at quantificationlib) - [TODO] add ensemble methods SC-MQ, MC-SQ, MC-MQ - [TODO] add HistNetQ - [TODO] add CDE-iteration and Bayes-CDE methods diff --git a/docs/build/html/_modules/quapy/method/aggregative.html b/docs/build/html/_modules/quapy/method/aggregative.html index c99b3be..09ac4fa 100644 --- a/docs/build/html/_modules/quapy/method/aggregative.html +++ b/docs/build/html/_modules/quapy/method/aggregative.html @@ -778,8 +778,13 @@ [docs] class CC(AggregativeCrispQuantifier): """ - The most basic Quantification method. One that simply classifies all instances and counts how many have been - attributed to each of the classes in order to compute class prevalence estimates. + `Classify & Count` (CC), the most basic quantification method, one that + simply classifies all instances and counts how many have been attributed to + each class in order to compute class prevalence estimates. This baseline is + the unadjusted estimator discussed, among others, in + `Forman, G. (2008). Quantifying counts and costs via classification. + Data Mining and Knowledge Discovery, 17, 164-206 + <https://link.springer.com/article/10.1007/s10618-008-0097-y>`_. :param classifier: a sklearn's Estimator that generates a classifier """ @@ -816,8 +821,13 @@ [docs] class PCC(AggregativeSoftQuantifier): """ - `Probabilistic Classify & Count <https://ieeexplore.ieee.org/abstract/document/5694031>`_, - the probabilistic variant of CC that relies on the posterior probabilities returned by a probabilistic classifier. + `Probabilistic Classify & Count` (PCC), the probabilistic variant of CC + that relies on the posterior probabilities returned by a probabilistic + classifier, introduced in + `Bella, A., Ferri, C., Hernández-Orallo, J., and Ramírez-Quintana, M.J. + (2010). Quantification via probability estimators. In Proceedings of the + 2010 IEEE International Conference on Data Mining (ICDM 2010) + <https://ieeexplore.ieee.org/abstract/document/5694031>`_. :param classifier: a sklearn's Estimator that generates a classifier """ @@ -849,9 +859,12 @@ [docs] class ACC(AggregativeCrispQuantifier): """ - `Adjusted Classify & Count <https://link.springer.com/article/10.1007/s10618-008-0097-y>`_, - the "adjusted" variant of :class:`CC`, that corrects the predictions of CC - according to the `misclassification rates`. + `Adjusted Classify & Count` (ACC), the "adjusted" variant of :class:`CC` + that corrects the predictions of CC according to the + misclassification rates, originally proposed in + `Forman, G. (2008). Quantifying counts and costs via classification. + Data Mining and Knowledge Discovery, 17, 164-206 + <https://link.springer.com/article/10.1007/s10618-008-0097-y>`_. :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']` @@ -1008,8 +1021,13 @@ [docs] class PACC(AggregativeSoftQuantifier): """ - `Probabilistic Adjusted Classify & Count <https://ieeexplore.ieee.org/abstract/document/5694031>`_, - the probabilistic variant of ACC that relies on the posterior probabilities returned by a probabilistic classifier. + `Probabilistic Adjusted Classify & Count` (PACC), the probabilistic + variant of ACC that relies on the posterior probabilities returned by a + probabilistic classifier, introduced in + `Bella, A., Ferri, C., Hernández-Orallo, J., and Ramírez-Quintana, M.J. + (2010). Quantification via probability estimators. In Proceedings of the + 2010 IEEE International Conference on Data Mining (ICDM 2010) + <https://ieeexplore.ieee.org/abstract/document/5694031>`_. :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']` @@ -2279,6 +2297,10 @@ KDEyHD = _kdey.KDEyHD KDEyCS = _kdey.KDEyCS +from . import _edy + +EDy = _edy.EDy + # --------------------------------------------------------------- # aliases # --------------------------------------------------------------- @@ -2290,6 +2312,7 @@ ExpectationMaximizationQuantifier = EMQ SLD = EMQ DistributionMatchingY = DMy +EnergyDistanceY = EDy HellingerDistanceY = HDy HistoricalHDy = DMy.HDy MedianSweep = MS diff --git a/docs/source/index.md b/docs/source/index.md index a9e4b64..3a78ba8 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -13,7 +13,7 @@ QuaPy is an open-source Python framework for quantification, also known as supervised prevalence estimation or learning to quantify. It is designed with research and experimental analysis in mind, and combines datasets, protocols, evaluation measures, visualization tools, and a broad collection of -quantification methods in a single coherent workflow. +quantification methods in a single workflow. ``` `````{grid} 1 1 2 2 diff --git a/docs/source/manuals/datasets.md b/docs/source/manuals/datasets.md index dff406c..50ed9ee 100644 --- a/docs/source/manuals/datasets.md +++ b/docs/source/manuals/datasets.md @@ -415,11 +415,12 @@ ECML-PKDD 2024, Vilnius, Lithuania. ## Image Embedding Datasets QuaPy also provides a collection of image datasets in the form of pre-generated -embeddings, hosted in [Zenodo](https://zenodo.org/records/21131944). +embeddings. These -embeddings were generated using [this extraction scripts](https://github.com/pglez82/visiondatasets_quapy). +embeddings were generated using [this extraction script](https://github.com/pglez82/visiondatasets_quapy) +and are hosted in [Zenodo](https://zenodo.org/records/21131944). -The current public interface is: +An example of current public interface is: ```python import quapy as qp @@ -432,21 +433,19 @@ data = qp.datasets.fetch_image_embeddings( train, test = data.train_test ``` -The available datasets are: +The available datasets are in `qp.datasets.IMAGE_DATASETS`, and include 6 datasets: -```python -qp.datasets.IMAGE_DATASETS -# ['cifar10', 'cifar100', 'cifar100coarse', 'svhn', 'fashionmnist', 'mnist'] -``` +* `cifar10`, `cifar100`, and `cifar100coarse`: + [Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. Technical report, University of Toronto, 2009.](https://cave.cs.toronto.edu/kriz/learning-features-2009-TR.pdf) +* `mnist`: + [Yann LeCun, Corinna Cortes, and Christopher J. C. Burges. The MNIST database of handwritten digits. 1998.](http://yann.lecun.com/exdb/mnist/) +* `fashionmnist`: + [Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv:1708.07747, 2017.](https://arxiv.org/abs/1708.07747) +* `svhn`: + [Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Baolin Wu, Andrew Y. Ng, et al. Reading digits in natural images with unsupervised feature learning. NIPS Workshop, 2011.](https://static.googleusercontent.com/media/research.google.com/es//pubs/archive/37648.pdf) -The available embedding types are: -```python -qp.datasets.IMAGE_EMBEDDINGS -# ['features', 'logits', 'predictions'] -``` - -where: +The available embedding types are in `qp.datasets.IMAGE_EMBEDDINGS`, and include: * `features` are the penultimate-layer representations * `logits` are the pre-activation outputs of the neural model @@ -456,49 +455,20 @@ The datasets correspond to frozen neural representations extracted from models trained on image classification tasks. QuaPy downloads them automatically on first use and stores them locally for fast reuse. -### Train/Test Semantics - Each dataset is internally organised into three splits: `train`, `val`, and `test`. The `train` split was used to train the neural model that produced the -embeddings, while `val` and `test` were not seen during neural training. - -For this reason, the default setting is: - -```python -data = qp.datasets.fetch_image_embeddings(..., heldout_only=True) -``` - -which returns: - -* `train = val` -* `test = test` - +embeddings, while `val` and `test` were not seen during neural training. +For this reason, the default setting indicates `heldout_only=True`, meaning +that the returned dataset will take the validation partition as the training +set, and the test partition as the test set. This is often the most convenient choice for quantification experiments, since it avoids training quantifiers on examples that were already used to train the embedding model. -If instead you want to use all the available non-test data, you can set: - -```python -data = qp.datasets.fetch_image_embeddings(..., heldout_only=False) -``` - -In this case, the returned training set is the union of the original neural +If instead you want to use all the available non-test data, you can set `heldout_only=False`, +in which case, the returned training set is the union of the original neural training split and the validation split. -### Sources - -The image datasets currently available through `fetch_image_embeddings` are: - -* cifar10, cifar100, cifar100coarse: - [Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. Technical report, University of Toronto, 2009.](https://cave.cs.toronto.edu/kriz/learning-features-2009-TR.pdf) -* mnist: - [Yann LeCun, Corinna Cortes, and Christopher J. C. Burges. The MNIST database of handwritten digits. 1998.](http://yann.lecun.com/exdb/mnist/) -* fashionmnist: - [Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv:1708.07747, 2017.](https://arxiv.org/abs/1708.07747) -* svhn: - [Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Baolin Wu, Andrew Y. Ng, et al. Reading digits in natural images with unsupervised feature learning. NIPS Workshop, 2011.](https://static.googleusercontent.com/media/research.google.com/es//pubs/archive/37648.pdf) - Some statistics are shown in the following table: | Dataset | backbone | classes | neural network train size | validation size | test size | feature dim | logit dim | prediction dim | type | diff --git a/docs/source/manuals/methods.md b/docs/source/manuals/methods.md index af5a54a..7851362 100644 --- a/docs/source/manuals/methods.md +++ b/docs/source/manuals/methods.md @@ -118,11 +118,15 @@ in evaluation. QuaPy implements the four CC variants, i.e.: * _CC_ (Classify & Count), the simplest aggregative quantifier; one that - simply relies on the label predictions of a classifier to deliver class estimates. -* _ACC_ (Adjusted Classify & Count), the adjusted variant of CC. + classifies all instances and computes the prevalence of the predicted labels. + This baseline is discussed, among others, in [Forman (2008)](https://link.springer.com/article/10.1007/s10618-008-0097-y). +* _ACC_ (Adjusted Classify & Count), the adjusted variant of CC, originally + proposed in [Forman (2008)](https://link.springer.com/article/10.1007/s10618-008-0097-y). * _PCC_ (Probabilistic Classify & Count), the probabilistic variant of CC that -relies on the soft estimations (or posterior probabilities) returned by a (probabilistic) classifier. -* _PACC_ (Probabilistic Adjusted Classify & Count), the adjusted variant of PCC. + relies on the posterior probabilities returned by a probabilistic classifier, + introduced in [Bella et al. (2010)](https://ieeexplore.ieee.org/abstract/document/5694031). +* _PACC_ (Probabilistic Adjusted Classify & Count), the adjusted variant of PCC, + also introduced in [Bella et al. (2010)](https://ieeexplore.ieee.org/abstract/document/5694031). The following code serves as a complete example using CC equipped with a SVM as the classifier: @@ -142,7 +146,7 @@ svm = LinearSVC() # (an alias is available in qp.method.aggregative.ClassifyAndCount) model = qp.method.aggregative.CC(svm) model.fit(Xtr, ytr) -estim_prevalence = model.predict(test.instances) +estim_prevalence = model.predict(test.X) ``` The same code could be used to instantiate an ACC, by simply replacing @@ -319,7 +323,21 @@ model.fit(*train.Xy) estim_prevalence = model.predict(test.X) ``` -### Hellinger Distance y (HDy) +### Distribution Matching + +Distribution Matching (DM) methods search for the mixture parameter (the sought class prevalence values) +yielding the mixture between the class-wise representations that best matches the test distribution. +Different criteria for deciding how this matching is assessed, and different ways for modelling the +distributions give rise to different instantiations of DM methods. + +The following methods are here discussed because they rely on a surrogate classifier for representing +the distributions, albeit different non-aggregative variants of them do often exist. Aside from this, +the formulation of DM methods is flexible enough as to accomodate methods that were proposed under a different +framework; examples include ACC and PACC. + +See the frameworks by [Firat](https://arxiv.org/abs/1606.00868), [Bunse](https://dl.gi.de/items/5a61f30f-6c84-4165-bd92-9098bd9e91aa), [Garg et al.](https://dl.acm.org/doi/10.5555/3495724.3496001), or [Dussap](https://theses.hal.science/tel-04931123), for more details. + +#### Hellinger Distance y (HDy) Implementation of the method based on the Hellinger Distance y (HDy) proposed by [González-Castro, V., Alaiz-Rodríguez, R., and Alegre, E. (2013). Class distribution @@ -354,23 +372,93 @@ model.fit(*dataset.training.Xy) estim_prevalence = model.predict(dataset.test.X) ``` -QuaPy also provides an implementation of the generalized -"Distribution Matching" approaches for multiclass, inspired by the framework -of [Firat (2016)](https://arxiv.org/abs/1606.00868). One can instantiate -a variant of HDy for multiclass quantification as follows: +#### Generalized Distribution Matching y (DMy) + +QuaPy also provides a generalized posterior-space distribution-matching +quantifier for binary or multiclass problems, implemented as +`qp.method.aggregative.DMy`. This class follows the generic distribution +matching view discussed by [Firat (2016)](https://arxiv.org/abs/1606.00868): +it represents class-conditional posterior distributions by histograms and then +searches for the prevalence vector whose mixture best matches the test +distribution. + +`DMy` is intentionally flexible and exposes three main design choices: the +number of histogram bins (`nbins`), the divergence to minimize (`divergence`, +e.g., `'HD'` or `'topsoe'`), and whether to match PDFs or CDFs (`cdf`). The +optimization routine can also be selected through `search`; the default +`'optim_minimize'` works for multiclass problems, while `'linear_search'` and +`'ternary_search'` are binary-only. A multiclass HDy-like instance can be +obtained as: ```python -mutliclassHDy = qp.method.aggregative.DMy(classifier=LogisticRegression(), divergence='HD', cdf=False) -``` +multiclass_hdy = qp.method.aggregative.DMy( + classifier=LogisticRegression(), + divergence='HD', + cdf=False, +) +``` -QuaPy also provides an implementation of the "DyS" -framework proposed by [Maletzke et al (2020)](https://ojs.aaai.org/index.php/AAAI/article/view/4376) -and the "SMM" method proposed by [Hassan et al (2019)](https://ieeexplore.ieee.org/document/9260028) -(thanks to _Pablo González_ for the contributions!) +#### DyS -A Bayesian distribution-matching counterpart is also available; see the +QuaPy implements the binary `DyS` framework proposed by +[Maletzke et al. (2020)](https://ojs.aaai.org/index.php/AAAI/article/view/4376) +as `qp.method.aggregative.DyS`. Conceptually, `DyS` can be seen as a +generalization of HDy in which the prevalence is found by ternary search over a +distribution-matching objective. In QuaPy, the user can select the number of +histogram bins (`n_bins`), the divergence (`divergence`), and the optimization +tolerance (`tol`). + +#### Energy Distance y (EDy) + +QuaPy also adapts `EDy` from [quantificationlib](https://github.com/AICGijon/quantificationlib), +which is available as `qp.method.aggregative.EDy`. + +This +method replaces histogram matching with an energy-distance formulation defined +directly on posterior-probability vectors and solves the resulting optimization +problem by quadratic programming. The method is proposed in +[Castaño et al.'s (2024)](https://ieeexplore.ieee.org/document/9791435/) paper. + +In QuaPy, `EDy` works for binary and +multiclass problems and lets the user choose the pairwise distance through the +`distance` parameter (`'manhattan'`, `'euclidean'`, or a custom callable). +Because the optimization relies on `quadprog`, this method requires the +optional dependency `pip install quadprog`. + +#### SMM + +QuaPy also includes the binary `SMM` method of +[Hassan et al. (2019)](https://ieeexplore.ieee.org/document/9260028), +available as `qp.method.aggregative.SMM`. This is a very lightweight +distribution-matching variant in which the posterior representation is reduced +to class-wise means rather than full histograms, making it conceptually close +to PACC. + + +#### Kernel Density Estimation methods (KDEy) + +QuaPy provides implementations for the three variants +of KDE-based methods proposed in +_[Moreo, A., González, P. and del Coz, J.J.. +Kernel Density Estimation for Multiclass Quantification. +Machine Learning. Vol 114 (92), 2025](https://link.springer.com/article/10.1007/s10994-024-06726-5)_ +(a [preprint](https://arxiv.org/abs/2401.00490) is available online). +The variants differ in the divergence metric to be minimized: + +- KDEy-HD: minimizes the (squared) Hellinger Distance and solves the problem via a Monte Carlo approach +- KDEy-CS: minimizes the Cauchy-Schwarz divergence and solves the problem via a closed-form solution +- KDEy-ML: minimizes the Kullback-Leibler divergence and solves the problem via maximum-likelihood + +These methods are specifically devised for multiclass problems (although they can tackle +binary problems too). + +All KDE-based methods depend on the hyperparameter `bandwidth` of the kernel. Typical values +that can be explored in model selection range in [0.01, 0.25]. Previous experiments reveal the methods' performance +varies smoothly at small variations of this hyperparameter. + +A Bayesian counterpart is available as well; see the {ref}`Bayesian Quantification Methods section ` -for `PQ` (Precise Quantifier). +for `BayesianKDEy`. ### Explicit Loss Minimization @@ -444,44 +532,21 @@ import quapy as qp from quapy.method.aggregative import SVMQ # load a single-label dataset (this one contains 3 classes) -dataset = qp.datasets.fetch_twitter('hcr', pickle=True) +train, test = qp.datasets.fetch_twitter('hcr', pickle=True).train_test # let qp know where svmperf is qp.environ['SVMPERF_HOME'] = '../svm_perf_quantification' model = newOneVsAll(SVMQ(), n_jobs=-1) # run them on parallel -model.fit(dataset.training) -estim_prevalence = model.predict(dataset.test.instances) +model.fit(*train.Xy) +estim_prevalence = model.predict(test.X) ``` Check the examples on [explicit loss minimization](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/17.explicit_loss_minimization.py) and on [one versus all quantification](https://github.com/HLT-ISTI/QuaPy/blob/devel/examples/10.one_vs_all.py) for more details. **Note** that the _one versus all_ approach is considered inappropriate under prior probability shift, though. -### Kernel Density Estimation methods (KDEy) -QuaPy provides implementations for the three variants -of KDE-based methods proposed in -_[Moreo, A., González, P. and del Coz, J.J.. -Kernel Density Estimation for Multiclass Quantification. -Machine Learning. Vol 114 (92), 2025](https://link.springer.com/article/10.1007/s10994-024-06726-5)_ -(a [preprint](https://arxiv.org/abs/2401.00490) is available online). -The variants differ in the divergence metric to be minimized: - -- KDEy-HD: minimizes the (squared) Hellinger Distance and solves the problem via a Monte Carlo approach -- KDEy-CS: minimizes the Cauchy-Schwarz divergence and solves the problem via a closed-form solution -- KDEy-ML: minimizes the Kullback-Leibler divergence and solves the problem via maximum-likelihood - -These methods are specifically devised for multiclass problems (although they can tackle -binary problems too). - -All KDE-based methods depend on the hyperparameter `bandwidth` of the kernel. Typical values -that can be explored in model selection range in [0.01, 0.25]. Previous experiments reveal the methods' performance -varies smoothly at small variations of this hyperparameter. - -A Bayesian counterpart is available as well; see the -{ref}`Bayesian Quantification Methods section ` -for `BayesianKDEy`. ## Non-Aggregative Methods diff --git a/docs/source/manuals/plots/app-grid.png b/docs/source/manuals/plots/app-grid.png new file mode 100644 index 0000000..d8cce7a Binary files /dev/null and b/docs/source/manuals/plots/app-grid.png differ diff --git a/docs/source/manuals/plots/app-kraemer.png b/docs/source/manuals/plots/app-kraemer.png new file mode 100644 index 0000000..1c842d1 Binary files /dev/null and b/docs/source/manuals/plots/app-kraemer.png differ diff --git a/docs/source/manuals/plots/dirichlet.png b/docs/source/manuals/plots/dirichlet.png new file mode 100644 index 0000000..664ffc6 Binary files /dev/null and b/docs/source/manuals/plots/dirichlet.png differ diff --git a/docs/source/manuals/plots/npp.png b/docs/source/manuals/plots/npp.png new file mode 100644 index 0000000..b43d382 Binary files /dev/null and b/docs/source/manuals/plots/npp.png differ diff --git a/docs/source/manuals/plots/simplex_visualization.png b/docs/source/manuals/plots/simplex_visualization.png new file mode 100644 index 0000000..2f6c570 Binary files /dev/null and b/docs/source/manuals/plots/simplex_visualization.png differ diff --git a/docs/source/manuals/plotting.md b/docs/source/manuals/plotting.md index cf685fe..153ccb4 100644 --- a/docs/source/manuals/plotting.md +++ b/docs/source/manuals/plotting.md @@ -284,4 +284,8 @@ qp.plot.plot_simplex( See the dedicated [example](https://github.com/HLT-ISTI/QuaPy/blob/master/examples/19.visualizing_simplex.py) -for a slightly richer illustration. +for a slightly richer illustration. The current example combines a posterior +cloud, the true/training/predicted prevalences, a smooth density surface, and +a region induced by Bonferroni-corrected 95% confidence intervals. + +![simplex visualization](./plots/simplex_visualization.png) diff --git a/docs/source/manuals/protocols.md b/docs/source/manuals/protocols.md index 17bc41a..aae89cb 100644 --- a/docs/source/manuals/protocols.md +++ b/docs/source/manuals/protocols.md @@ -43,21 +43,21 @@ desired prevalence values covering the full spectrum. In APP, the user specifies the number of (equally distant) points to be generated from the interval [0,1]; -in QuaPy this is achieved by setting _n_prevpoints_. -For example, if _n_prevpoints=11_ then, for each class, the prevalence values +in QuaPy this is achieved by setting _n_prevalences_. +For example, if _n_prevalences=11_ then, for each class, the prevalence values [0., 0.1, 0.2, ..., 1.] will be used. This means that, for two classes, the number of different prevalence values will be 11 (since, once the prevalence of one class is determined, the other one is constrained). For 3 classes, the number of valid combinations can be obtained as 11 + 10 + ... + 1 = 66. In general, the number of valid combinations that will be produced for a given -value of n_prevpoints can be consulted by invoking +value of _n_prevalences_ can be consulted by invoking _num_prevalence_combinations_, e.g.: ```python import quapy.functional as F -n_prevpoints = 21 +n_prevalences = 21 n_classes = 4 -n = F.num_prevalence_combinations(n_prevpoints, n_classes, n_repeats=1) +n = F.num_prevalence_combinations(n_prevalences, n_classes, n_repeats=1) ``` in this example, _n=1771_. Note the last argument, _n_repeats_, that @@ -74,13 +74,13 @@ _get_nprevpoints_approximation_, e.g.: ```python budget = 5000 -n_prevpoints = F.get_nprevpoints_approximation(budget, n_classes, n_repeats=1) -n = F.num_prevalence_combinations(n_prevpoints, n_classes, n_repeats=1) -print(f'by setting n_prevpoints={n_prevpoints} the number of evaluations for {n_classes} classes will be {n}') +n_prevalences = F.get_nprevpoints_approximation(budget, n_classes, n_repeats=1) +n = F.num_prevalence_combinations(n_prevalences, n_classes, n_repeats=1) +print(f'by setting n_prevalences={n_prevalences} the number of evaluations for {n_classes} classes will be {n}') ``` this will produce the following output: ``` -by setting n_prevpoints=30 the number of evaluations for 4 classes will be 4960 +by setting n_prevalences=30 the number of evaluations for 4 classes will be 4960 ``` The following code shows an example of usage of APP for model selection @@ -129,6 +129,15 @@ in such cases QuaPy takes the value of _qp.environ['SAMPLE_SIZE']_. This protocol is useful for testing a quantifier under conditions of _prior probability shift_. +The following ternary plot, generated by [example 21](https://github.com/HLT-ISTI/QuaPy/blob/master/examples/21.visualizing_protocols.py), +shows the prevalence values covered by a grid-based APP in a three-class problem (`academic-success`): + +![APP grid protocol](./plots/app-grid.png) + +Each point corresponds to one sampled prevalence vector. As expected, the +points lie on a regular grid over the simplex, ensuring systematic coverage of +the prevalence space. + ## Sampling from the unit-simplex, the Uniform-Prevalence Protocol (UPP) Generating all possible combinations from a grid of prevalence values (APP) in @@ -148,7 +157,7 @@ for sampling from the unit-simplex as many vectors of prevalence values as indic in the _repeats_ parameter. UPP can be instantiated as: ```python -protocol = qp.in_protocol.UPP(test, repeats=100) +protocol = qp.protocol.UPP(test, repeats=100) ``` This is the most convenient protocol for datasets @@ -157,6 +166,16 @@ containing many classes; see, e.g., and is useful for testing a quantifier under conditions of _prior probability shift_. +The next plot shows one such protocol, labelled in example 21 as +_APP(Kraemer)_, to emphasize that it plays the role of an artificial-prevalence +protocol without relying on a fixed grid: + +![UPP / APP Kraemer protocol](./plots/app-kraemer.png) + +Unlike grid-based APP, UPP does not force prevalence vectors to lie on a +regular lattice. Instead, it spreads samples over the simplex in a +statistically uniform way, making it attractive when the number of classes is +large and exhaustive grids become impractical. ## Natural-Prevalence Protocol @@ -168,9 +187,47 @@ All other things being equal, this protocol can be used just like APP or UPP, and is instantiated via: ```python -protocol = qp.in_protocol.NPP(test, repeats=100) +protocol = qp.protocol.NPP(test, repeats=100) ``` +The prevalence coverage of NPP is much more concentrated, since the samples are +obtained by plain random subsampling from the test set and therefore remain +close to its natural prevalence: + +![Natural-prevalence protocol](./plots/npp.png) + +This makes NPP useful when one wants to evaluate quantifiers under mild or +realistic drift conditions, but much less suitable than APP or UPP for stress +testing performance across the full simplex. + +## Dirichlet Protocol + +QuaPy also implements a :class:`DirichletProtocol`, which samples prevalence +vectors from a Dirichlet distribution before drawing the corresponding sample +from the labelled collection: + +```python +protocol = qp.protocol.DirichletProtocol(test, alpha=0.2, repeats=100) +``` + +The parameter `alpha` controls how concentrated the protocol is. Small values +of `alpha` favour sparse prevalence vectors near the corners of the simplex, +while larger values generate more balanced mixtures. When all entries of +`alpha` are equal to 1, the protocol becomes uniformly distributed over the +simplex, similarly in spirit to UPP. + +The following plot shows the effect of a sparse prior with `alpha=0.2`: + +![Dirichlet protocol](./plots/dirichlet.png) + +Compared to UPP, the mass is clearly pulled towards the vertices and edges, +thus producing more extreme label-shift scenarios. + +The parameter `alpha` in the Dirichlet distribution is typically defined as an +array of shape `(n_classes)`. When the user specifies a single value, QuaPy +broadcasts this value for all classes. Conversely, a different value can be +specified for each class. + ## Other protocols Other protocols exist in QuaPy and will be added to the `qp.protocol.py` module. \ No newline at end of file diff --git a/examples/21.visualizing_protocols.py b/examples/21.visualizing_protocols.py new file mode 100644 index 0000000..8f87e69 --- /dev/null +++ b/examples/21.visualizing_protocols.py @@ -0,0 +1,51 @@ +import numpy as np +import quapy as qp +from quapy.data.datasets import fetch_UCIMulticlassDataset +from quapy.protocol import APP, NPP, UPP, DirichletProtocol + + +""" +Ternary plots showcasing different sampling protocols. +""" + +rng = np.random.default_rng(0) + +train, test = fetch_UCIMulticlassDataset(dataset_name='academic-success').train_test + +train_prev = { + 'points': train.prevalence(), + 'label': 'training prevalence', + 'style': {'s': 70, 'color': 'darkorange'}, +} + +def protocols(): + yield 'app-grid', 'Artificial Prevalence Protocol (grid)', APP(test, n_prevalences=21, repeats=1, sample_size=100) + yield 'app-kraemer', 'Artificial Prevalence Protocol (Kraemer)', UPP(test, repeats=5000, sample_size=500) + yield 'npp', 'Natural Prevalence Protocol', NPP(test, repeats=1000, sample_size=100) + yield 'dirichlet', 'Dirichlet(alpha=0.2)', DirichletProtocol(test, alpha=0.2, repeats=5000, sample_size=100) + +for file_name, prot_name, protocol in protocols(): + app_points = { + 'points': [prev for _, prev in protocol()], + 'label': prot_name, + 'style': {'s': 15, 'alpha': 0.5, 'color': 'steelblue', 'edgecolors': 'none'}, + } + + point_layers = [ + app_points, + train_prev, + ] + + dispersion = 0.1 + + qp.plot.plot_simplex( + point_layers=point_layers, + class_names=['class A', 'class B', 'class C'], + #title='Ternary prevalence visualisation', + legend_ncol=3, + figsize=(7.2, 5.8), + class_name_fontsize=9, + title_fontsize=10, + legend_fontsize=8, + savepath=f'./plots/{file_name}.png', + ) diff --git a/quapy/method/__init__.py b/quapy/method/__init__.py index 7f05027..1273a5a 100644 --- a/quapy/method/__init__.py +++ b/quapy/method/__init__.py @@ -26,6 +26,7 @@ AGGREGATIVE_METHODS = { aggregative.MS, aggregative.MS2, aggregative.DMy, + aggregative.EDy, aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, @@ -55,6 +56,7 @@ MULTICLASS_METHODS = { aggregative.PACC, aggregative.RLLS, aggregative.EMQ, + aggregative.EDy, aggregative.KDEyML, aggregative.KDEyCS, aggregative.KDEyHD, diff --git a/quapy/method/_edy.py b/quapy/method/_edy.py new file mode 100644 index 0000000..9d0f938 --- /dev/null +++ b/quapy/method/_edy.py @@ -0,0 +1,269 @@ +from typing import Callable, Union + +import numpy as np +from sklearn.base import BaseEstimator +from sklearn.metrics.pairwise import euclidean_distances, manhattan_distances + +import quapy as qp +import quapy.functional as F +from quapy.method._helper import _get_quadprog +from quapy.method.aggregative import AggregativeSoftQuantifier + + +class EDy(AggregativeSoftQuantifier): + """ + Energy Distance y (EDy), a posterior-space distribution-matching quantifier + based on energy distance. + + The method represents each class by the posterior-probability vectors + produced by a probabilistic classifier on validation data, and estimates the + test prevalence vector by matching the test posterior distribution against + the class-conditional validation distributions through an energy-distance + objective solved as a quadratic program. The method is therefore another + instance of the general mixture-matching view of quantification, but it + operates directly on posterior vectors rather than on histogram summaries. + + This implementation works for binary and multiclass single-label + quantification and relies on the optional ``quadprog`` dependency. It was + adapted to QuaPy's current aggregative API from the original implementation + available in `quantificationlib `_. + + The current implementation follows the energy-distance formulation discussed + in: + + * Alberto Castaño, Laura Morán-Fernández, Jaime Alonso, + Verónica Bolón-Canedo, Amparo Alonso-Betanzos, and Juan José del Coz. + *An analysis of quantification methods based on matching distributions*. + * Hideko Kawakubo, Marthinus Christoffel du Plessis, and Masashi Sugiyama + (2016). *Computationally efficient class-prior estimation under class + balance change using energy distance*. IEICE Transactions on Information + and Systems, 99(1):176-186. + + :param classifier: a scikit-learn ``BaseEstimator``, or ``None`` to use + ``qp.environ['DEFAULT_CLS']`` + :param fit_classifier: whether to train the learner (default ``True``). + Set to ``False`` if the learner has already been trained outside the + quantifier + :param val_split: specification of the data used for generating validation + posterior probabilities. This can be an integer (default ``5``) for + k-fold cross-validation, a float in ``(0, 1)`` for a held-out split, + or a tuple ``(X, y)`` with explicit validation data + :param distance: distance used to compare posterior vectors. Valid string + aliases are ``'manhattan'`` (default) and ``'euclidean'``; a custom + callable compatible with pairwise-distance signatures can also be used + :param n_jobs: number of parallel workers (default ``None``, meaning the + value is taken from the environment) + """ + + def __init__( + self, + classifier: BaseEstimator = None, + fit_classifier: bool = True, + val_split=5, + distance: Union[str, Callable] = 'manhattan', + n_jobs=None, + ): + super().__init__(classifier, fit_classifier, val_split) + self.distance = distance + self.n_jobs = qp._get_njobs(n_jobs) + self.train_n_cls_i_ = None + self.train_distrib_ = None + self.K_ = None + self.G_ = None + self.C_ = None + self.b_ = None + self.a_ = None + + def _check_init_parameters(self): + self.distance = self._resolve_distance_function(self.distance) + + @staticmethod + def _resolve_distance_function(distance): + if isinstance(distance, str): + if distance == 'manhattan': + return manhattan_distances + if distance == 'euclidean': + return euclidean_distances + raise ValueError( + f"unknown distance {distance!r}; valid aliases are 'manhattan' and 'euclidean'" + ) + if not hasattr(distance, '__call__'): + raise ValueError('distance must be a valid string alias or a callable function') + return distance + + def _is_pd(self, m): + """Check whether a symmetric matrix is positive definite. + + This helper is used before invoking ``quadprog`` because the quadratic + term of the optimization problem must be positive definite. + """ + return self._dpofa(m)[0] == 0 + + def _dpofa(self, m): + """Factor a symmetric positive definite matrix. + + This is a lightweight Python adaptation of the ``dpofa`` routine used by + ``quadprog``. Here it is mainly employed as a numerical check while + preparing the quadratic-program matrix. + """ + r = np.array(m, copy=True) + n = len(r) + for k in range(n): + s = 0.0 + if k >= 1: + for i in range(k): + t = r[i, k] + if i > 0: + t = t - np.sum(r[0:i, i] * r[0:i, k]) + t = t / r[i, i] + r[i, k] = t + s = s + t * t + s = r[k, k] - s + if s <= 0.0: + return k + 1, r + r[k, k] = np.sqrt(s) + return 0, r + + def _nearest_pd(self, A): + """Project a matrix onto the cone of positive-definite matrices. + + In some cases the matrix induced by the energy-distance objective is not + numerically positive definite, even though the underlying optimization + problem is well posed. In those cases we replace it with the nearest + positive-definite approximation before calling ``quadprog``. + """ + B = (A + A.T) / 2 + _, s, V = np.linalg.svd(B) + H = V.T @ np.diag(s) @ V + A2 = (B + H) / 2 + A3 = (A2 + A2.T) / 2 + + if self._is_pd(A3): + return A3 + + spacing = np.spacing(np.linalg.norm(A)) + identity_matrix = np.eye(A.shape[0]) + k = 1 + while not self._is_pd(A3): + mineig = np.min(np.real(np.linalg.eigvals(A3))) + A3 += identity_matrix * (-mineig * k ** 2 + spacing) + k += 1 + + return A3 + + def _compute_ed_param_train(self, distance_func, train_distrib, classes, n_cls_i): + """Pre-compute the training-side terms of the ED optimization problem. + + Given the class-conditional posterior clouds observed on validation + data, this routine computes the pairwise average distances between + classes and derives the matrices required by the quadratic program. + These terms depend only on the validation distribution and can therefore + be cached after ``aggregation_fit``. + """ + n_classes = len(classes) + K = np.zeros((n_classes, n_classes), dtype=float) + for i in range(n_classes): + K[i, i] = distance_func(train_distrib[classes[i]], train_distrib[classes[i]]).sum() + for j in range(i + 1, n_classes): + K[i, j] = distance_func(train_distrib[classes[i]], train_distrib[classes[j]]).sum() + K[j, i] = K[i, j] + + K = K / np.dot(n_cls_i, n_cls_i.T) + + B = np.zeros((n_classes - 1, n_classes - 1), dtype=float) + for i in range(n_classes - 1): + B[i, i] = -K[i, i] - K[-1, -1] + 2 * K[i, -1] + for j in range(n_classes - 1): + if j == i: + continue + B[i, j] = -K[i, j] - K[-1, -1] + K[i, -1] + K[j, -1] + + G = 2 * B + if not self._is_pd(G): + G = self._nearest_pd(G) + + C = -np.vstack([np.ones((1, n_classes - 1)), -np.eye(n_classes - 1)]).T + b = -np.array([1] + [0] * (n_classes - 1), dtype=float) + + return K, G, C, b + + def _compute_ed_param_test(self, distance_func, train_distrib, test_distrib, K, classes, n_cls_i): + """Compute the test-dependent linear term of the ED objective. + + Once the training-side matrices have been computed, each new test sample + only requires estimating the distances between its posterior cloud and + the class-conditional validation clouds. + """ + n_classes = len(classes) + Kt = np.zeros(n_classes, dtype=float) + for i in range(n_classes): + Kt[i] = distance_func(train_distrib[classes[i]], test_distrib).sum() + + Kt = Kt / (n_cls_i.squeeze() * float(len(test_distrib))) + return 2 * (-Kt[:-1] + K[:-1, -1] + Kt[-1] - K[-1, -1]) + + def _solve_ed(self, G, a, C, b): + """Solve the energy-distance quadratic program. + + The optimization is carried out over the first ``n_classes - 1`` + prevalences; the prevalence of the last class is recovered afterwards by + the simplex constraint. The resulting vector is finally normalized as a + precaution against small numerical deviations. + """ + quadprog = _get_quadprog() + sol = quadprog.solve_qp(G=G, a=a, C=C, b=b) + prevalences = sol[0] + prevalences = np.append(prevalences, 1 - prevalences.sum()) + return F.normalize_prevalence(prevalences, method='clip') + + def aggregation_fit(self, classif_predictions, labels): + """ + Estimate the class-conditional posterior distributions on validation + data and pre-compute the quadratic-program parameters that depend only + on the training side. + + In EDy, the validation posteriors are not discretized into histograms. + Instead, each class is represented by the cloud of posterior vectors + observed for that class, and these clouds are then compared through the + selected pairwise distance. + + :param classif_predictions: posterior probabilities returned by the + classifier on validation data + :param labels: true labels associated to each posterior vector + """ + posteriors = np.asarray(classif_predictions, dtype=float) + labels = np.asarray(labels) + + self.train_distrib_ = { + class_: posteriors[labels == class_] for class_ in self.classes_ + } + self.train_n_cls_i_ = np.asarray( + [[len(self.train_distrib_[class_])] for class_ in self.classes_], + dtype=float, + ) + + self.K_, self.G_, self.C_, self.b_ = self._compute_ed_param_train( + self.distance, + self.train_distrib_, + self.classes_, + self.train_n_cls_i_, + ) + return self + + def aggregate(self, posteriors: np.ndarray): + """Estimate the prevalence vector for a test sample. + + :param posteriors: posterior probabilities returned by the classifier + for the instances in the test sample + :return: a prevalence vector of shape ``(n_classes,)`` + """ + posteriors = np.asarray(posteriors, dtype=float) + self.a_ = self._compute_ed_param_test( + self.distance, + self.train_distrib_, + posteriors, + self.K_, + self.classes_, + self.train_n_cls_i_, + ) + return self._solve_ed(G=self.G_, a=self.a_, C=self.C_, b=self.b_) diff --git a/quapy/method/_helper.py b/quapy/method/_helper.py index 0e5792a..48d0cb6 100644 --- a/quapy/method/_helper.py +++ b/quapy/method/_helper.py @@ -32,6 +32,16 @@ def _get_cvxpy(): return cp +def _get_quadprog(): + try: + import quadprog + except ImportError as exc: + raise ImportError( + "EDy requires the optional 'quadprog' package." + ) from exc + return quadprog + + def _labels_to_indices(labels, classes): encoder = LabelEncoder().fit(classes) return encoder.transform(labels) diff --git a/quapy/method/aggregative.py b/quapy/method/aggregative.py index 5e57ee4..2aa160c 100644 --- a/quapy/method/aggregative.py +++ b/quapy/method/aggregative.py @@ -371,8 +371,13 @@ class BinaryAggregativeQuantifier(AggregativeQuantifier, BinaryQuantifier): # ------------------------------------ class CC(AggregativeCrispQuantifier): """ - The most basic Quantification method. One that simply classifies all instances and counts how many have been - attributed to each of the classes in order to compute class prevalence estimates. + `Classify & Count` (CC), the most basic quantification method, one that + simply classifies all instances and counts how many have been attributed to + each class in order to compute class prevalence estimates. This baseline is + the unadjusted estimator discussed, among others, in + `Forman, G. (2008). Quantifying counts and costs via classification. + Data Mining and Knowledge Discovery, 17, 164-206 + `_. :param classifier: a sklearn's Estimator that generates a classifier """ @@ -400,8 +405,13 @@ class CC(AggregativeCrispQuantifier): class PCC(AggregativeSoftQuantifier): """ - `Probabilistic Classify & Count `_, - the probabilistic variant of CC that relies on the posterior probabilities returned by a probabilistic classifier. + `Probabilistic Classify & Count` (PCC), the probabilistic variant of CC + that relies on the posterior probabilities returned by a probabilistic + classifier, introduced in + `Bella, A., Ferri, C., Hernández-Orallo, J., and Ramírez-Quintana, M.J. + (2010). Quantification via probability estimators. In Proceedings of the + 2010 IEEE International Conference on Data Mining (ICDM 2010) + `_. :param classifier: a sklearn's Estimator that generates a classifier """ @@ -424,9 +434,12 @@ class PCC(AggregativeSoftQuantifier): class ACC(AggregativeCrispQuantifier): """ - `Adjusted Classify & Count `_, - the "adjusted" variant of :class:`CC`, that corrects the predictions of CC - according to the `misclassification rates`. + `Adjusted Classify & Count` (ACC), the "adjusted" variant of :class:`CC` + that corrects the predictions of CC according to the + misclassification rates, originally proposed in + `Forman, G. (2008). Quantifying counts and costs via classification. + Data Mining and Knowledge Discovery, 17, 164-206 + `_. :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']` @@ -568,8 +581,13 @@ class ACC(AggregativeCrispQuantifier): class PACC(AggregativeSoftQuantifier): """ - `Probabilistic Adjusted Classify & Count `_, - the probabilistic variant of ACC that relies on the posterior probabilities returned by a probabilistic classifier. + `Probabilistic Adjusted Classify & Count` (PACC), the probabilistic + variant of ACC that relies on the posterior probabilities returned by a + probabilistic classifier, introduced in + `Bella, A., Ferri, C., Hernández-Orallo, J., and Ramírez-Quintana, M.J. + (2010). Quantification via probability estimators. In Proceedings of the + 2010 IEEE International Conference on Data Mining (ICDM 2010) + `_. :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']` @@ -1712,6 +1730,10 @@ KDEyML = _kdey.KDEyML KDEyHD = _kdey.KDEyHD KDEyCS = _kdey.KDEyCS +from . import _edy + +EDy = _edy.EDy + # --------------------------------------------------------------- # aliases # --------------------------------------------------------------- @@ -1723,6 +1745,7 @@ ProbabilisticAdjustedClassifyAndCount = PACC ExpectationMaximizationQuantifier = EMQ SLD = EMQ DistributionMatchingY = DMy +EnergyDistanceY = EDy HellingerDistanceY = HDy HistoricalHDy = DMy.HDy MedianSweep = MS diff --git a/quapy/tests/test_methods.py b/quapy/tests/test_methods.py index 54954a0..bb9b060 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, HDx -from quapy.method.aggregative import ACC, DMy, KDEyCS, RLLS +from quapy.method.aggregative import ACC, DMy, EDy, KDEyCS, RLLS from quapy.method.meta import Ensemble from quapy.functional import check_prevalence_vector from quapy.tests._synthetic import make_dataset @@ -20,6 +20,7 @@ OPTIONAL_AGGREGATIVE_METHODS = { 'BayesianMAPLS', 'PQ', 'RLLS', + 'EDy', } @@ -139,6 +140,19 @@ class TestMethods(unittest.TestCase): self.assertTrue(check_prevalence_vector(estim_prevalences)) + def test_edy(self): + try: + import quadprog # noqa: F401 + except ImportError: + return + + dataset = TestMethods.tiny_dataset_multiclass + q = EDy(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_dmy_noncanonical_labels(self): dataset = TestMethods.tiny_dataset_multiclass label_names = np.asarray(['class-a', 'class-c', 'class-z'])