integrating bayesian methods and related functionality, plus unit test refactor

CHANGE_LOG.txt

@@ -2,8 +2,15 @@ Change Log 0.2.1
 -----------------
 
 - Improved documentation of confidence regions.
+- Added Bayesian KDEy and Bayesian MAPLS quantifiers.
+- Added temperature calibration utilities for Bayesian confidence-aware methods.
+- Added compositional CLR and ILR transformations.
+- Extended KDEy with Aitchison/ILR kernels, shrinkage, and improved numerical stability.
+- Added TemperatureScalingFromLogits for calibrating pretrained logits.
+- Added DirichletProtocol for prevalence sampling from Dirichlet priors.
 - Added ReadMe method by Daniel Hopkins and Gary King
 - Internal index in LabelledCollection is now "lazy", and is only constructed if required.
+- Improved unit testing and separated integration tests
 
 Change Log 0.2.0
 -----------------
@@ -214,4 +221,3 @@ Change Log 0.1.7
     any instance of BaseQuantifier), and a subclass of it called OneVsAllAggregative which implements the
     classify / aggregate interface. Both are instances of OneVsAll. There is a method getOneVsAll that returns the
     best instance based on the type of quantifier.
-

quapy/__init__.py

@@ -7,12 +7,16 @@ from . import functional
 from . import method
 from . import evaluation
 from . import protocol
-from . import plot
 from . import util
 from . import model_selection
 from . import classification
 import os
 
+try:
+    from . import plot
+except ImportError:
+    plot = None
+
 __version__ = '0.2.1'
 
 
@@ -74,4 +78,3 @@ def _get_classifier(classifier):
         raise ValueError('neither classifier nor qp.environ["DEFAULT_CLS"] have been specified')
     return classifier
 
-

quapy/classification/__init__.py

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

quapy/classification/calibration.py

@@ -1,8 +1,9 @@
 from copy import deepcopy
 
-from abstention.calibration import NoBiasVectorScaling, TempScaling, VectorScaling
 from sklearn.base import BaseEstimator, clone
 from sklearn.model_selection import cross_val_predict, train_test_split
+from sklearn.preprocessing import LabelEncoder
+from sklearn.utils.validation import check_X_y
 import numpy as np
 
 
@@ -11,6 +12,17 @@ import numpy as np
 # see https://github.com/kundajelab/abstention
 
 
+def _require_abstention_calibration():
+    try:
+        from abstention.calibration import NoBiasVectorScaling, TempScaling, VectorScaling
+    except ImportError as exc:
+        raise ImportError(
+            "Calibration methods in quapy.classification.calibration require the optional "
+            "'abstention' package."
+        ) from exc
+    return NoBiasVectorScaling, TempScaling, VectorScaling
+
+
 class RecalibratedProbabilisticClassifier:
     """
     Abstract class for (re)calibration method from `abstention.calibration`, as defined in
@@ -142,6 +154,7 @@ class NBVSCalibration(RecalibratedProbabilisticClassifierBase):
     """
 
     def __init__(self, classifier, val_split=5, n_jobs=None, verbose=False):
+        NoBiasVectorScaling, _, _ = _require_abstention_calibration()
         self.classifier = classifier
         self.calibrator = NoBiasVectorScaling(verbose=verbose)
         self.val_split = val_split
@@ -164,6 +177,7 @@ class BCTSCalibration(RecalibratedProbabilisticClassifierBase):
     """
 
     def __init__(self, classifier, val_split=5, n_jobs=None, verbose=False):
+        _, TempScaling, _ = _require_abstention_calibration()
         self.classifier = classifier
         self.calibrator = TempScaling(verbose=verbose, bias_positions='all')
         self.val_split = val_split
@@ -186,6 +200,7 @@ class TSCalibration(RecalibratedProbabilisticClassifierBase):
     """
 
     def __init__(self, classifier, val_split=5, n_jobs=None, verbose=False):
+        _, TempScaling, _ = _require_abstention_calibration()
         self.classifier = classifier
         self.calibrator = TempScaling(verbose=verbose)
         self.val_split = val_split
@@ -208,9 +223,84 @@ class VSCalibration(RecalibratedProbabilisticClassifierBase):
     """
 
     def __init__(self, classifier, val_split=5, n_jobs=None, verbose=False):
+        _, _, VectorScaling = _require_abstention_calibration()
         self.classifier = classifier
         self.calibrator = VectorScaling(verbose=verbose)
         self.val_split = val_split
         self.n_jobs = n_jobs
         self.verbose = verbose
 
+
+class TemperatureScalingFromLogits(BaseEstimator):
+    """
+    Calibrates a matrix of logits by learning a temperature-scaling mapping
+    with the calibration methods from `abstention.calibration`.
+
+    This estimator is useful when the inputs are already logits produced by a
+    pretrained classifier, and the goal is to transform them directly into
+    calibrated posterior probabilities without retraining the underlying model.
+
+    :param bias_corrected: if True, uses Bias-Corrected Temperature Scaling
+        (BCTS); otherwise, uses standard Temperature Scaling (TS)
+    :param verbose: whether the underlying calibrator should display progress
+        information
+    """
+
+    def __init__(self, bias_corrected=False, verbose=False):
+        self.bias_corrected = bias_corrected
+        self.verbose = verbose
+
+    def fit(self, X, y):
+        """
+        Fits the logits calibrator.
+
+        :param X: array-like of shape `(n_samples, n_classes)` containing
+            logits
+        :param y: array-like of shape `(n_samples,)` containing class labels
+        :return: self
+        """
+        X, y = check_X_y(X, y)
+
+        self.label_encoder_ = LabelEncoder()
+        y_enc = self.label_encoder_.fit_transform(y)
+        self.classes_ = self.label_encoder_.classes_
+
+        n_classes = len(self.classes_)
+        logits_dim = X.shape[1]
+        if n_classes != logits_dim:
+            raise ValueError(
+                f'mismatch between the number of classes ({n_classes}) and the '
+                f'dimensionality of the logits ({logits_dim})'
+            )
+
+        _, TempScaling, _ = _require_abstention_calibration()
+        calibrator = TempScaling(
+            verbose=self.verbose,
+            bias_positions='all' if self.bias_corrected else [],
+        )
+        self.calibrator_ = calibrator
+        self.calibration_function_ = calibrator(X, np.eye(n_classes)[y_enc])
+        return self
+
+    def predict_proba(self, X):
+        """
+        Converts logits into calibrated posterior probabilities.
+
+        :param X: array-like of shape `(n_samples, n_classes)` containing
+            logits
+        :return: array-like of shape `(n_samples, n_classes)` with calibrated
+            posterior probabilities
+        """
+        return self.calibration_function_(X)
+
+    def predict(self, X):
+        """
+        Predicts class labels after calibration.
+
+        :param X: array-like of shape `(n_samples, n_classes)` containing
+            logits
+        :return: array-like of shape `(n_samples,)` with class label
+            predictions
+        """
+        posteriors = self.predict_proba(X)
+        return self.label_encoder_.inverse_transform(np.argmax(posteriors, axis=1))

quapy/classification/methods.py

@@ -1,3 +1,4 @@
+import numpy as np
 from sklearn.base import BaseEstimator
 from sklearn.decomposition import TruncatedSVD
 from sklearn.linear_model import LogisticRegression
@@ -95,3 +96,23 @@ class LowRankLogisticRegression(BaseEstimator):
         if self.pca is None:
             return X
         return self.pca.transform(X)
+
+
+class MockClassifierFromPosteriors(BaseEstimator):
+    """
+    Mock classifier that bypasses classifier training when the input instances
+    are already posterior probabilities produced by a pretrained probabilistic
+    classifier.
+
+    :param X: arrays of shape `(n_samples, n_classes)` are interpreted as posterior probabilities
+    """
+
+    def fit(self, X, y):
+        self.classes_ = np.sort(np.unique(y))
+        return self
+
+    def predict(self, X):
+        return np.argmax(X, axis=1)
+
+    def predict_proba(self, X):
+        return X

quapy/data/datasets.py

@@ -3,7 +3,6 @@ from contextlib import contextmanager
 import zipfile
 from os.path import join
 import pandas as pd
-from ucimlrepo import fetch_ucirepo
 from quapy.data.base import Dataset, LabelledCollection
 from quapy.data.preprocessing import text2tfidf, reduce_columns
 from quapy.data.preprocessing import standardize as standardizer
@@ -12,6 +11,17 @@ from quapy.util import download_file_if_not_exists, download_file, get_quapy_hom
 from sklearn.preprocessing import StandardScaler
 
 
+def _fetch_ucirepo(*args, **kwargs):
+    try:
+        from ucimlrepo import fetch_ucirepo
+    except ImportError as exc:
+        raise ImportError(
+            "UCI dataset fetching requires the optional 'ucimlrepo' package. "
+            "Install it to use fetch_UCIBinaryDataset or fetch_UCIMulticlassDataset."
+        ) from exc
+    return fetch_ucirepo(*args, **kwargs)
+
+
 REVIEWS_SENTIMENT_DATASETS = ['hp', 'kindle', 'imdb']
 
 TWITTER_SENTIMENT_DATASETS_TEST = [
@@ -486,7 +496,7 @@ def fetch_UCIBinaryLabelledCollection(dataset_name, data_home=None, standardize=
         # fall back to direct download when needed
         if group == "german":
             with download_tmp_file("statlog/german", "german.data-numeric") as tmp:
-                df = pd.read_csv(tmp, header=None, sep="\\s+")
+                df = pd.read_csv(tmp, header=None, delim_whitespace=True)
             X, y = df.iloc[:, 0:24].astype(float).values, df[24].astype(int).values
         elif group == "ctg":
             with download_tmp_file("00193", "CTG.xls") as tmp:
@@ -500,11 +510,11 @@ def fetch_UCIBinaryLabelledCollection(dataset_name, data_home=None, standardize=
             y = df["NSP"].astype(int).values
         elif group == "semeion":
             with download_tmp_file("semeion", "semeion.data") as tmp:
-                df = pd.read_csv(tmp, header=None, sep="\\s+")
+                df = pd.read_csv(tmp, header=None, sep='\\s+')
             X = df.iloc[:, 0:256].astype(float).values
             y = df[263].values  # 263 stands for digit 8 (labels are one-hot vectors from col 256-266)
         else:
-            df = fetch_ucirepo(id=id)
+            df = _fetch_ucirepo(id=id)
             X, y = df.data.features.to_numpy(), df.data.targets.to_numpy().squeeze()
 
         # transform data when needed before returning (returned data will be pickled)
@@ -616,8 +626,8 @@ def fetch_UCIMulticlassDataset(
         are taken for training, and the rest (irrespective of `min_test_split`) is taken for test.
     :param max_train_instances: maximum number of instances to keep for training (defaults to 25000);
         set to -1 or None to avoid this check
-    :param min_class_support: minimum number of istances per class. Classes with fewer instances
-        are discarded (deafult is 100)
+    :param min_class_support: integer or float, the minimum number or proportion of istances per class.
+        Classes with fewer instances are discarded (deafult is 100).
     :param standardize: indicates whether the covariates should be standardized or not (default is True). If requested,
         standardization applies after the LabelledCollection is split, that is, the mean an std are computed only on the
         training portion of the data.
@@ -673,6 +683,11 @@ def fetch_UCIMulticlassLabelledCollection(dataset_name, data_home=None, min_clas
         f'Name {dataset_name} does not match any known dataset from the ' \
         f'UCI Machine Learning datasets repository (multiclass). ' \
         f'Valid ones are {UCI_MULTICLASS_DATASETS}'
+
+    assert (min_class_support is None or
+            ((isinstance(min_class_support, int) and min_class_support >= 0) or
+             (isinstance(min_class_support, float) and 0. <= min_class_support < 1.))), \
+        f'invalid value for {min_class_support=}; expected non negative integer or float in [0,1)'
     
     if data_home is None:
         data_home = get_quapy_home()
@@ -739,24 +754,41 @@ def fetch_UCIMulticlassLabelledCollection(dataset_name, data_home=None, min_clas
 
     file = join(data_home, 'uci_multiclass', dataset_name+'.pkl')
     
+    def dummify_categorical_features(df_features, dataset_id):
+        categorical_features = {
+            158: ["S1", "C1", "S2", "C2", "S3", "C3", "S4", "C4", "S5", "C5"],  # poker_hand
+        }
+
+        categorical = categorical_features.get(dataset_id, [])
+
+        X = df_features.copy()
+        if categorical:
+            X[categorical] = X[categorical].astype("category")
+            X = pd.get_dummies(X, columns=categorical, drop_first=True)
+
+        return X
+
     def download(id, name):
-        df = fetch_ucirepo(id=id)
+        df = _fetch_ucirepo(id=id)
 
-        df.data.features = pd.get_dummies(df.data.features, drop_first=True)
-        X, y = df.data.features.to_numpy(dtype=np.float64), df.data.targets.to_numpy().squeeze()
+        X_df = dummify_categorical_features(df.data.features, id)
+        X = X_df.to_numpy(dtype=np.float64)
+        y = df.data.targets.to_numpy().squeeze()
 
-        assert y.ndim == 1, 'more than one y'
+        assert y.ndim == 1, f'error: the dataset {id=} {name=} has more than one target variable'
 
         classes = np.sort(np.unique(y))
         y = np.searchsorted(classes, y)
         return LabelledCollection(X, y)
 
-    def filter_classes(data: LabelledCollection, min_ipc):
-        if min_ipc is None:
-            min_ipc = 0
+    def filter_classes(data: LabelledCollection, min_class_support):
+        if min_class_support is None or min_class_support == 0.:
+            return data
+        if isinstance(min_class_support, float):
+            min_class_support = int(len(data) * min_class_support)
         classes = data.classes_
-        # restrict classes to only those with at least min_ipc instances
-        classes = classes[data.counts() >= min_ipc]
+        # restrict classes to only those with at least min_class_support instances
+        classes = classes[data.counts() >= min_class_support]
         # filter X and y keeping only datapoints belonging to valid classes
         filter_idx = np.isin(data.y, classes)
         X, y = data.X[filter_idx], data.y[filter_idx]

quapy/error.py

@@ -128,6 +128,78 @@ def se(prevs_true, prevs_hat):
     return ((prevs_hat - prevs_true) ** 2).mean(axis=-1)
 
 
+def sre(prevs_true, prevs_hat, prevs_train, eps=0.):
+    """
+    Computes the squared ratio error between two prevalence vectors.
+    The squared ratio error between prevalence vectors :math:`p` and
+    :math:`\\hat{p}` with training prevalence :math:`p^{tr}` is:
+    :math:`SRE(p,\\hat{p},p^{tr})=\\frac{1}{|\\mathcal{Y}|}\\sum_{i \\in \\mathcal{Y}}(w_i-\\hat{w}_i)^2`,
+    where :math:`w_i=\\frac{p_i}{p^{tr}_i}`.
+
+    :param prevs_true: array-like with the true prevalence values
+    :param prevs_hat: array-like with the predicted prevalence values
+    :param prevs_train: array-like with the training prevalence values, or a single
+        prevalence vector when all comparisons refer to the same training set
+    :param eps: smoothing factor for the prevalence values (default 0, i.e., no smoothing)
+    :return: squared ratio error
+    """
+    prevs_true = np.asarray(prevs_true)
+    prevs_hat = np.asarray(prevs_hat)
+    prevs_train = np.asarray(prevs_train)
+    assert prevs_true.shape == prevs_hat.shape, f'wrong shape {prevs_true.shape=} vs {prevs_hat.shape=}'
+    assert prevs_true.shape[-1] == prevs_train.shape[-1], 'wrong shape for training prevalence'
+    if prevs_true.ndim == 2 and prevs_train.ndim == 1:
+        prevs_train = np.tile(prevs_train, reps=(prevs_true.shape[0], 1))
+    if eps > 0:
+        prevs_true = smooth(prevs_true, eps)
+        prevs_hat = smooth(prevs_hat, eps)
+        prevs_train = smooth(prevs_train, eps)
+
+    n_classes = prevs_true.shape[-1]
+    w = prevs_true / prevs_train
+    w_hat = prevs_hat / prevs_train
+    return (1. / n_classes) * np.sum((w - w_hat) ** 2., axis=-1)
+
+
+def msre(prevs_true, prevs_hat, prevs_train, eps=0.):
+    """
+    Computes the mean squared ratio error (see :meth:`quapy.error.sre`) across the sample pairs.
+
+    :param prevs_true: array-like of shape `(n_samples, n_classes,)` with the true prevalence values
+    :param prevs_hat: array-like of shape equal to prevs_true with the predicted prevalence values
+    :param prevs_train: array-like with the training prevalence values
+    :param eps: smoothing factor (default 0, i.e., no smoothing)
+    :return: mean squared ratio error
+    """
+    return np.mean(sre(prevs_true, prevs_hat, prevs_train, eps))
+
+
+def aitchisondist(prevs_true, prevs_hat):
+    """
+    Computes the Aitchison distance between two prevalence vectors.
+
+    :param prevs_true: array-like with the true prevalence values
+    :param prevs_hat: array-like with the predicted prevalence values
+    :return: Aitchison distance
+    """
+    from quapy.functional import CLRtransformation
+
+    clr = CLRtransformation()
+    return np.linalg.norm(clr(prevs_true) - clr(prevs_hat), axis=-1)
+
+
+def maitchisondist(prevs_true, prevs_hat):
+    """
+    Computes the mean Aitchison distance (see :meth:`quapy.error.aitchisondist`)
+    across the sample pairs.
+
+    :param prevs_true: array-like with the true prevalence values
+    :param prevs_hat: array-like with the predicted prevalence values
+    :return: mean Aitchison distance
+    """
+    return np.mean(aitchisondist(prevs_true, prevs_hat))
+
+
 def mkld(prevs_true, prevs_hat, eps=None):
     """Computes the mean Kullback-Leibler divergence (see :meth:`quapy.error.kld`) across the
     sample pairs. The distributions are smoothed using the `eps` factor
@@ -374,8 +446,8 @@ def __check_eps(eps=None):
 
 
 CLASSIFICATION_ERROR = {f1e, acce}
-QUANTIFICATION_ERROR = {mae, mnae, mrae, mnrae, mse, mkld, mnkld}
-QUANTIFICATION_ERROR_SINGLE = {ae, nae, rae, nrae, se, kld, nkld}
+QUANTIFICATION_ERROR = {mae, mnae, mrae, mnrae, mse, mkld, mnkld, msre, maitchisondist}
+QUANTIFICATION_ERROR_SINGLE = {ae, nae, rae, nrae, se, kld, nkld, sre, aitchisondist}
 QUANTIFICATION_ERROR_SMOOTH = {kld, nkld, rae, nrae, mkld, mnkld, mrae}
 CLASSIFICATION_ERROR_NAMES = {func.__name__ for func in CLASSIFICATION_ERROR}
 QUANTIFICATION_ERROR_NAMES = {func.__name__ for func in QUANTIFICATION_ERROR}
@@ -387,6 +459,9 @@ ERROR_NAMES = \
 f1_error = f1e
 acc_error = acce
 mean_absolute_error = mae
+squared_ratio_error = sre
+dist_aitchison = aitchisondist
+mean_dist_aitchison = maitchisondist
 absolute_error = ae
 mean_relative_absolute_error = mrae
 relative_absolute_error = rae

quapy/functional.py

@@ -1,11 +1,15 @@
 import warnings
+from abc import ABC, abstractmethod
 from collections import defaultdict
+from functools import lru_cache
 from typing import Literal, Union, Callable
 from numpy.typing import ArrayLike
 
 import scipy
 import numpy as np
 
+import quapy as qp
+
 
 # ------------------------------------------------------------------------------------------
 # General utils
@@ -649,3 +653,105 @@ def solve_adjustment(
         raise ValueError(f'unknown {solver=}')
 
 
+# ------------------------------------------------------------------------------------------
+# Transformations from Compositional analysis
+# ------------------------------------------------------------------------------------------
+
+class CompositionalTransformation(ABC):
+    """
+    Abstract class of transformations for compositional data.
+    """
+
+    EPSILON = 1e-12
+
+    @abstractmethod
+    def __call__(self, X):
+        ...
+
+    @abstractmethod
+    def inverse(self, Z):
+        ...
+
+
+class CLRtransformation(CompositionalTransformation):
+    """
+    Centered log-ratio (CLR) transformation.
+    """
+
+    def __call__(self, X):
+        X = np.asarray(X)
+        X = qp.error.smooth(X, self.EPSILON)
+        geometric_mean = np.exp(np.mean(np.log(X), axis=-1, keepdims=True))
+        return np.log(X / geometric_mean)
+
+    def inverse(self, Z):
+        return scipy.special.softmax(Z, axis=-1)
+
+
+class ILRtransformation(CompositionalTransformation):
+    """
+    Isometric log-ratio (ILR) transformation.
+    """
+
+    def __call__(self, X):
+        X = np.asarray(X)
+        X = qp.error.smooth(X, self.EPSILON)
+        basis = self.get_V(X.shape[-1])
+        return np.log(X) @ basis.T
+
+    def inverse(self, Z):
+        Z = np.asarray(Z)
+        basis = self.get_V(Z.shape[-1] + 1)
+        logp = Z @ basis
+        p = np.exp(logp)
+        return p / np.sum(p, axis=-1, keepdims=True)
+
+    @lru_cache(maxsize=None)
+    def get_V(self, k):
+        helmert = np.zeros((k, k))
+        for i in range(1, k):
+            helmert[i, :i] = 1
+            helmert[i, i] = -i
+            helmert[i] = helmert[i] / np.sqrt(i * (i + 1))
+        return helmert[1:, :]
+
+
+def normalized_entropy(p):
+    """
+    Computes the normalized Shannon entropy of a prevalence vector.
+
+    :param p: array-like prevalence vector summing to 1
+    :return: float in [0,1]
+    """
+    p = np.asarray(p)
+    entropy = scipy.stats.entropy(p)
+    max_entropy = np.log(len(p))
+    return np.clip(entropy / max_entropy, 0, 1)
+
+
+def antagonistic_prevalence(p, strength=1):
+    """
+    Reflects a prevalence vector in ILR space and maps it back to the simplex.
+
+    :param p: array-like prevalence vector
+    :param strength: reflection strength in ILR space
+    :return: prevalence vector in the simplex
+    """
+    ilr = ILRtransformation()
+    z = ilr(p)
+    z_ant = -strength * z
+    return ilr.inverse(z_ant)
+
+
+def in_simplex(x, atol=1e-8):
+    """
+    Checks whether points lie in the probability simplex.
+
+    :param x: array-like of shape `(n_classes,)` or `(n_points, n_classes)`
+    :param atol: numerical tolerance for the unit-sum check
+    :return: boolean or boolean array
+    """
+    x = np.asarray(x)
+    non_negative = np.all(x >= 0, axis=-1)
+    sum_to_one = np.isclose(x.sum(axis=-1), 1.0, atol=atol)
+    return non_negative & sum_to_one

quapy/method/__init__.py

@@ -3,6 +3,7 @@ from sklearn.exceptions import ConvergenceWarning
 warnings.simplefilter("ignore", ConvergenceWarning)
 
 from . import confidence
+from . import _bayesian
 from . import base
 from . import aggregative
 from . import non_aggregative
@@ -29,6 +30,8 @@ AGGREGATIVE_METHODS = {
     aggregative.KDEyHD,
     # aggregative.OneVsAllAggregative,
     confidence.BayesianCC,
+    _bayesian.BayesianKDEy,
+    _bayesian.BayesianMAPLS,
     confidence.PQ,
 }
 
@@ -53,7 +56,9 @@ MULTICLASS_METHODS = {
     aggregative.KDEyML,
     aggregative.KDEyCS,
     aggregative.KDEyHD,
-    confidence.BayesianCC
+    confidence.BayesianCC,
+    _bayesian.BayesianKDEy,
+    _bayesian.BayesianMAPLS,
 }
 
 NON_AGGREGATIVE_METHODS = {
@@ -71,4 +76,3 @@ QUANTIFICATION_METHODS = AGGREGATIVE_METHODS | NON_AGGREGATIVE_METHODS | META_ME
 
 
 
-

quapy/method/_bayesian.py

@@ -1,22 +1,49 @@
 """
-Utility functions for `Bayesian quantification <https://arxiv.org/abs/2302.09159>`_ methods.
+Utilities and methods for Bayesian quantification.
 """
-import numpy as np
+import contextlib
+import copy
 import importlib.resources
+import logging
+import os
+import sys
+from collections.abc import Iterable
+from numbers import Number, Real
+
+import numpy as np
+from joblib import Parallel, delayed
+from sklearn.base import BaseEstimator
+from tqdm import tqdm
+
+import quapy as qp
+import quapy.functional as F
+from quapy.data import LabelledCollection
+from quapy.method._kdey import KDEBase
+from quapy.method.aggregative import AggregativeSoftQuantifier
+from quapy.method.confidence import ConfidenceRegionABC, WithConfidenceABC
+from quapy.protocol import AbstractProtocol
 
 try:
     import jax
     import jax.numpy as jnp
+    import jax.random as jrandom
+    from jax.scipy.special import logsumexp as jax_logsumexp
     import numpyro
     import numpyro.distributions as dist
+    from numpyro.infer import MCMC, NUTS
     import stan
+    import stan.common
 
     DEPENDENCIES_INSTALLED = True
 except ImportError:
     jax = None
     jnp = None
+    jrandom = None
+    jax_logsumexp = None
     numpyro = None
     dist = None
+    MCMC = None
+    NUTS = None
     stan = None
 
     DEPENDENCIES_INSTALLED = False
@@ -27,28 +54,124 @@ P_TEST_C: str = "P_test(C)"
 P_C_COND_Y: str = "P(C|Y)"
 
 
-def model(n_c_unlabeled: np.ndarray, n_y_and_c_labeled: np.ndarray) -> None:
-    """
-    Defines a probabilistic model in `NumPyro <https://num.pyro.ai/>`_.
+def _require_bayesian_dependencies():
+    if not DEPENDENCIES_INSTALLED:
+        raise ImportError(
+            "Auxiliary dependencies are required. "
+            "Run `$ pip install quapy[bayes]` to install them."
+        )
 
-    :param n_c_unlabeled: a `np.ndarray` of shape `(n_predicted_classes,)`
-        with entry `c` being the number of instances predicted as class `c`.
-    :param n_y_and_c_labeled: a `np.ndarray` of shape `(n_classes, n_predicted_classes)`
-        with entry `(y, c)` being the number of instances labeled as class `y` and predicted as class `c`.
+
+def _resolve_dirichlet_prior(prior, n_classes, *, allow_mapls_priors=False, n_test=None, map_prev=None, map_lambda=None):
+    if isinstance(prior, str):
+        if prior == 'uniform':
+            return np.ones(n_classes, dtype=float)
+        if allow_mapls_priors and prior in {'map', 'map2'}:
+            if n_test is None or map_prev is None:
+                raise ValueError('MAPLS priors require n_test and map_prev')
+            if prior == 'map':
+                lam = map_lambda
+            else:
+                lam = get_lambda(
+                    test_probs=map_prev["test_probs"],
+                    pz=map_prev["train_prev"],
+                    q_prior=map_prev["map_estimate"],
+                    dvg=kl_div,
+                )
+            alpha_0 = alpha0_from_lambda(lam, n_test=n_test, n_classes=n_classes)
+            return np.full(n_classes, alpha_0, dtype=float)
+        raise ValueError(f"unknown prior specification {prior!r}")
+    if isinstance(prior, Number):
+        return np.full(n_classes, float(prior), dtype=float)
+
+    alpha = np.asarray(prior, dtype=float)
+    if alpha.ndim != 1 or len(alpha) != n_classes:
+        raise ValueError(f'wrong shape for prior; expected {n_classes} values, found shape {alpha.shape}')
+    return alpha
+
+
+def _validate_temperature(temperature):
+    if not isinstance(temperature, Real) or temperature <= 0:
+        raise ValueError(f'expected a positive real value for temperature; found {temperature!r}')
+    return float(temperature)
+
+
+def model_bayesianCC(
+    n_c_unlabeled: np.ndarray,
+    n_y_and_c_labeled: np.ndarray,
+    temperature: float,
+    alpha: np.ndarray,
+) -> None:
+    """
+    NumPyro model for BayesianCC.
     """
     n_y_labeled = n_y_and_c_labeled.sum(axis=1)
 
-    K = len(n_c_unlabeled)
-    L = len(n_y_labeled)
+    n_pred_classes = len(n_c_unlabeled)
+    n_classes = len(n_y_labeled)
 
-    pi_ = numpyro.sample(P_TEST_Y, dist.Dirichlet(jnp.ones(L)))
-    p_c_cond_y = numpyro.sample(P_C_COND_Y, dist.Dirichlet(jnp.ones(K).repeat(L).reshape(L, K)))
+    pi_ = numpyro.sample(P_TEST_Y, dist.Dirichlet(jnp.asarray(alpha, dtype=jnp.float32)))
+    p_c_cond_y = numpyro.sample(
+        P_C_COND_Y,
+        dist.Dirichlet(jnp.ones(n_pred_classes).repeat(n_classes).reshape(n_classes, n_pred_classes)),
+    )
 
-    with numpyro.plate('plate', L):
-        numpyro.sample('F_yc', dist.Multinomial(n_y_labeled, p_c_cond_y), obs=n_y_and_c_labeled)
+    if temperature == 1.0:
+        with numpyro.plate('plate', n_classes):
+            numpyro.sample('F_yc', dist.Multinomial(n_y_labeled, p_c_cond_y), obs=n_y_and_c_labeled)
+
+        p_c = numpyro.deterministic(P_TEST_C, jnp.einsum("yc,y->c", p_c_cond_y, pi_))
+        numpyro.sample('N_c', dist.Multinomial(jnp.sum(n_c_unlabeled), p_c), obs=n_c_unlabeled)
+        return
+
+    with numpyro.plate('plate_y', n_classes):
+        logp_f = dist.Multinomial(n_y_labeled, p_c_cond_y).log_prob(n_y_and_c_labeled)
+
+    numpyro.factor('F_yc_loglik', jnp.sum(logp_f) / temperature)
 
     p_c = numpyro.deterministic(P_TEST_C, jnp.einsum("yc,y->c", p_c_cond_y, pi_))
-    numpyro.sample('N_c', dist.Multinomial(jnp.sum(n_c_unlabeled), p_c), obs=n_c_unlabeled)
+    logp_n = dist.Multinomial(jnp.sum(n_c_unlabeled), p_c).log_prob(n_c_unlabeled)
+    numpyro.factor('N_c_loglik', logp_n / temperature)
+
+
+def model(n_c_unlabeled: np.ndarray, n_y_and_c_labeled: np.ndarray) -> None:
+    """
+    Backward-compatible BayesianCC model with a uniform prior.
+    """
+    alpha = np.ones(n_y_and_c_labeled.shape[0], dtype=float)
+    return model_bayesianCC(n_c_unlabeled, n_y_and_c_labeled, temperature=1.0, alpha=alpha)
+
+
+def sample_posterior_bayesianCC(
+    n_c_unlabeled: np.ndarray,
+    n_y_and_c_labeled: np.ndarray,
+    num_warmup: int,
+    num_samples: int,
+    alpha: np.ndarray,
+    temperature: float = 1.0,
+    seed: int = 0,
+) -> dict:
+    """
+    Samples from the BayesianCC posterior using NumPyro.
+    """
+    _require_bayesian_dependencies()
+    temperature = _validate_temperature(temperature)
+
+    mcmc = numpyro.infer.MCMC(
+        numpyro.infer.NUTS(model_bayesianCC),
+        num_warmup=num_warmup,
+        num_samples=num_samples,
+        progress_bar=False,
+    )
+    rng_key = jax.random.PRNGKey(seed)
+    mcmc.run(
+        rng_key,
+        n_c_unlabeled=n_c_unlabeled,
+        n_y_and_c_labeled=n_y_and_c_labeled,
+        temperature=temperature,
+        alpha=alpha,
+    )
+    return mcmc.get_samples()
 
 
 def sample_posterior(
@@ -59,77 +182,679 @@ def sample_posterior(
     seed: int = 0,
 ) -> dict:
     """
-    Samples from the Bayesian quantification model in NumPyro using the
-    `NUTS <https://arxiv.org/abs/1111.4246>`_ sampler.
-
-    :param n_c_unlabeled: a `np.ndarray` of shape `(n_predicted_classes,)`
-        with entry `c` being the number of instances predicted as class `c`.
-    :param n_y_and_c_labeled: a `np.ndarray` of shape `(n_classes, n_predicted_classes)`
-        with entry `(y, c)` being the number of instances labeled as class `y` and predicted as class `c`.
-    :param num_warmup: the number of warmup steps.
-    :param num_samples: the number of samples to draw.
-    :seed: the random seed.
-    :return: a `dict` with the samples. The keys are the names of the latent variables.
+    Backward-compatible wrapper around BayesianCC sampling.
     """
-    mcmc = numpyro.infer.MCMC(
-        numpyro.infer.NUTS(model),
+    alpha = np.ones(n_y_and_c_labeled.shape[0], dtype=float)
+    return sample_posterior_bayesianCC(
+        n_c_unlabeled=n_c_unlabeled,
+        n_y_and_c_labeled=n_y_and_c_labeled,
         num_warmup=num_warmup,
         num_samples=num_samples,
-        progress_bar=False
+        alpha=alpha,
+        temperature=1.0,
+        seed=seed,
     )
-    rng_key = jax.random.PRNGKey(seed)
-    mcmc.run(rng_key, n_c_unlabeled=n_c_unlabeled, n_y_and_c_labeled=n_y_and_c_labeled)
-    return mcmc.get_samples()
-
 
 
 def load_stan_file():
     return importlib.resources.files('quapy.method').joinpath('stan/pq.stan').read_text(encoding='utf-8')
 
+
+@contextlib.contextmanager
+def _suppress_stan_logging():
+    with open(os.devnull, "w") as devnull:
+        old_stderr = sys.stderr
+        sys.stderr = devnull
+        try:
+            yield
+        finally:
+            sys.stderr = old_stderr
+
+
 def pq_stan(stan_code, n_bins, pos_hist, neg_hist, test_hist, number_of_samples, num_warmup, stan_seed):
     """
-    Perform Bayesian prevalence estimation using a Stan model for probabilistic quantification.
-
-    This function builds and samples from a Stan model that implements a bin-based Bayesian
-    quantifier. It uses the class-conditional histograms of the classifier
-    outputs for positive and negative examples, along with the test histogram, to estimate
-    the posterior distribution of prevalence in the test set.
-
-    Parameters
-    ----------
-    stan_code : str
-        The Stan model code as a string. 
-    n_bins : int
-        Number of bins used to build the histograms for positive and negative examples.
-    pos_hist : array-like of shape (n_bins,)
-        Histogram counts of the classifier outputs for the positive class.
-    neg_hist : array-like of shape (n_bins,)
-        Histogram counts of the classifier outputs for the negative class.
-    test_hist : array-like of shape (n_bins,)
-        Histogram counts of the classifier outputs for the test set, binned using the same bins.
-    number_of_samples : int
-        Number of post-warmup samples to draw from the Stan posterior.
-    num_warmup : int
-        Number of warmup iterations for the sampler.
-    stan_seed : int
-        Random seed for Stan model compilation and sampling, ensuring reproducibility.
-
-    Returns
-    -------
-    prev_samples : numpy.ndarray
-        An array of posterior samples of the prevalence (`prev`) in the test set.
-        Each element corresponds to one draw from the posterior distribution.
+    Samples posterior prevalences for PQ from a Stan model.
     """
+    _require_bayesian_dependencies()
+    logging.getLogger("stan.common").setLevel(logging.ERROR)
 
     stan_data = {
-            'n_bucket': n_bins,
-            'train_neg': neg_hist.tolist(),
-            'train_pos': pos_hist.tolist(),
-            'test': test_hist.tolist(),
-            'posterior': 1
-        }
+        'n_bucket': n_bins,
+        'train_neg': neg_hist.tolist(),
+        'train_pos': pos_hist.tolist(),
+        'test': test_hist.tolist(),
+        'posterior': 1,
+    }
 
-    stan_model = stan.build(stan_code, data=stan_data, random_seed=stan_seed)
-    fit = stan_model.sample(num_chains=1, num_samples=number_of_samples,num_warmup=num_warmup)
+    with _suppress_stan_logging():
+        stan_model = stan.build(stan_code, data=stan_data, random_seed=stan_seed)
+        fit = stan_model.sample(num_chains=1, num_samples=number_of_samples, num_warmup=num_warmup)
 
     return fit['prev']
+
+
+class BayesianKDEy(AggregativeSoftQuantifier, KDEBase, WithConfidenceABC):
+    """
+    Bayesian version of KDEy.
+
+    This method relies on extra dependencies, which have to be installed via:
+    `$ pip install quapy[bayes]`
+
+    :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 classifier, or consider it
+        already fit
+    :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. Set to None when the method does not require any
+        validation data, in order to avoid that some portion of the training
+        data be wasted.
+    :param kernel: kernel function for KDE. Available kernels include
+        {'gaussian', 'aitchison', 'ilr'} (default 'gaussian')
+    :param bandwidth: bandwidth for the kernel (default 0.1)
+    :param shrinkage: regularization strength for Aitchison/ILR kernels
+        (default 0.0)
+    :param num_warmup: number of warmup iterations for the MCMC sampler
+        (default 500)
+    :param num_samples: number of samples to draw from the posterior
+        (default 1000)
+    :param mcmc_seed: random seed for the MCMC sampler (default 0)
+    :param confidence_level: float in [0,1] to construct a confidence region
+        around the point estimate (default 0.95)
+    :param region: string, set to `intervals` for constructing confidence
+        intervals (default), or to `ellipse` for constructing an ellipse in
+        the probability simplex, or to `ellipse-clr` for constructing an
+        ellipse in the Centered-Log Ratio (CLR) unconstrained space.
+    :param temperature: temperature (>0) for posterior calibration
+        (default 1.)
+    :param prior: an array-like with the alpha parameters of a Dirichlet
+        prior, a scalar real value to be broadcast to all classes, or the
+        string 'uniform' for a uniform, uninformative prior (default)
+    :param verbose: bool, whether to display the progress bar
+    """
+
+    def __init__(
+        self,
+        classifier: BaseEstimator = None,
+        fit_classifier=True,
+        val_split: int = 5,
+        kernel='gaussian',
+        bandwidth=0.1,
+        shrinkage=0.0,
+        num_warmup: int = 500,
+        num_samples: int = 1_000,
+        mcmc_seed: int = 0,
+        confidence_level: float = 0.95,
+        region: str = 'intervals',
+        temperature: float = 1.0,
+        prior='uniform',
+        verbose: bool = False,
+    ):
+        _require_bayesian_dependencies()
+        if num_warmup <= 0:
+            raise ValueError(f'parameter {num_warmup=} must be a positive integer')
+        if num_samples <= 0:
+            raise ValueError(f'parameter {num_samples=} must be a positive integer')
+
+        self.kernel = KDEBase._check_kernel(kernel)
+        self.bandwidth = KDEBase._check_bandwidth(bandwidth, self.kernel)
+        assert 0 <= shrinkage < 1, 'shrinkage must be in [0,1)'
+        assert self.kernel != 'gaussian' or shrinkage == 0, \
+            'shrinkage is only supported for Aitchison/ILR kernels'
+
+        super().__init__(classifier, fit_classifier, val_split)
+        self.shrinkage = float(shrinkage)
+        self.num_warmup = num_warmup
+        self.num_samples = num_samples
+        self.mcmc_seed = mcmc_seed
+        self.confidence_level = confidence_level
+        self.region = region
+        self.temperature = _validate_temperature(temperature)
+        self.prior = prior
+        self.verbose = verbose
+        self.prevalence_samples = None
+
+    def aggregation_fit(self, classif_predictions, labels):
+        self.mix_densities = self.get_mixture_components(
+            classif_predictions,
+            labels,
+            self.classes_,
+            self.bandwidth,
+            self.kernel,
+        )
+        return self
+
+    def sample_from_posterior(self, classif_predictions):
+        test_log_densities = np.asarray(
+            [self.pdf(kde_i, classif_predictions, self.kernel, log_densities=True) for kde_i in self.mix_densities]
+        )
+        alpha = _resolve_dirichlet_prior(self.prior, len(self.mix_densities))
+
+        mcmc = MCMC(
+            NUTS(self._numpyro_model),
+            num_warmup=self.num_warmup,
+            num_samples=self.num_samples,
+            num_chains=1,
+            progress_bar=self.verbose,
+        )
+        mcmc.run(jrandom.PRNGKey(self.mcmc_seed), test_log_densities=test_log_densities, alpha=alpha)
+        self.prevalence_samples = np.asarray(mcmc.get_samples()["prev"])
+        return self.prevalence_samples
+
+    def aggregate(self, classif_predictions: np.ndarray):
+        return self.sample_from_posterior(classif_predictions).mean(axis=0)
+
+    def predict_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
+        confidence_level = self.confidence_level if confidence_level is None else confidence_level
+        classif_predictions = self.classify(instances)
+        point_estimate = self.aggregate(classif_predictions)
+        region = WithConfidenceABC.construct_region(
+            self.prevalence_samples,
+            confidence_level=confidence_level,
+            method=self.region,
+        )
+        return point_estimate, region
+
+    def _numpyro_model(self, test_log_densities, alpha):
+        prev = numpyro.sample("prev", dist.Dirichlet(jnp.asarray(alpha)))
+        log_likelihood = jnp.sum(jax_logsumexp(jnp.log(prev)[:, None] + test_log_densities, axis=0))
+        numpyro.factor("loglik", (1.0 / self.temperature) * log_likelihood)
+
+
+class _JaxILRTransformation(F.CompositionalTransformation):
+    """
+    JAX-backed ILR transform used inside Bayesian MAPLS.
+    """
+
+    def __call__(self, X):
+        X = jnp.asarray(X)
+        X = qp.error.smooth(np.asarray(X), self.EPSILON)
+        X = jnp.asarray(X)
+        basis = jnp.asarray(self.get_V(X.shape[-1]))
+        return jnp.log(X) @ basis.T
+
+    def inverse(self, Z):
+        Z = jnp.asarray(Z)
+        basis = jnp.asarray(self.get_V(Z.shape[-1] + 1))
+        logp = Z @ basis
+        p = jnp.exp(logp)
+        return p / jnp.sum(p, axis=-1, keepdims=True)
+
+    def get_V(self, k):
+        return F.ILRtransformation().get_V(k)
+
+
+class BayesianMAPLS(AggregativeSoftQuantifier, WithConfidenceABC):
+    """
+    Bayesian variant of the MLLS/EMQ method proposed by
+    Ye, Changkun, et al. "Label shift estimation for class-imbalance problem:
+    A bayesian approach." Proceedings of the IEEE/CVF Winter Conference on
+    Applications of Computer Vision. 2024.
+
+    Code adapted from:
+    https://github.com/ChangkunYe/MAPLS/blob/main/label_shift/mapls.py
+
+    This method relies on extra dependencies, which have to be installed via:
+    `$ pip install quapy[bayes]`
+
+    :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 classifier, or consider it
+        already fit
+    :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. Set to None when the method does not require any
+        validation data, in order to avoid that some portion of the training
+        data be wasted.
+    :param exact_train_prev: set to True (default) for using the true training
+        prevalence as the initial observation; set to False for computing the
+        training prevalence as an estimate of it, i.e., as the expected value
+        of the posterior probabilities of the training instances.
+    :param num_warmup: number of warmup iterations for the MCMC sampler
+        (default 500)
+    :param num_samples: number of samples to draw from the posterior
+        (default 1000)
+    :param mcmc_seed: random seed for the MCMC sampler (default 0)
+    :param confidence_level: float in [0,1] to construct a confidence region
+        around the point estimate (default 0.95)
+    :param region: string, set to `intervals` for constructing confidence
+        intervals (default), or to `ellipse` for constructing an ellipse in
+        the probability simplex, or to `ellipse-clr` for constructing an
+        ellipse in the Centered-Log Ratio (CLR) unconstrained space.
+    :param temperature: temperature (>0) for posterior calibration
+        (default 1.)
+    :param prior: an array-like with the alpha parameters of a Dirichlet
+        prior, a scalar real value to be broadcast to all classes, or one of
+        {'uniform', 'map', 'map2'} (default 'uniform')
+    :param mapls_chain_init: whether to initialize the Markov chain with a
+        preliminary EM point estimate (default True)
+    :param verbose: bool, whether to display the progress bar
+        (default False)
+    """
+
+    def __init__(
+        self,
+        classifier: BaseEstimator = None,
+        fit_classifier=True,
+        val_split: int = 5,
+        exact_train_prev=True,
+        num_warmup: int = 500,
+        num_samples: int = 1_000,
+        mcmc_seed: int = 0,
+        confidence_level: float = 0.95,
+        region: str = 'intervals',
+        temperature: float = 1.0,
+        prior='uniform',
+        mapls_chain_init=True,
+        verbose=False,
+    ):
+        _require_bayesian_dependencies()
+        if num_warmup <= 0:
+            raise ValueError(f'parameter {num_warmup=} must be a positive integer')
+        if num_samples <= 0:
+            raise ValueError(f'parameter {num_samples=} must be a positive integer')
+        if not (
+            (isinstance(prior, str) and prior in {'uniform', 'map', 'map2'})
+            or isinstance(prior, Number)
+            or (isinstance(prior, Iterable) and all(isinstance(v, Number) for v in prior))
+        ):
+            raise ValueError(
+                f'wrong type for {prior=}; expected one of {{"uniform", "map", "map2"}}, '
+                'a real scalar, or an array-like of real values'
+            )
+
+        super().__init__(classifier, fit_classifier, val_split)
+        self.exact_train_prev = exact_train_prev
+        self.num_warmup = num_warmup
+        self.num_samples = num_samples
+        self.mcmc_seed = mcmc_seed
+        self.confidence_level = confidence_level
+        self.region = region
+        self.temperature = _validate_temperature(temperature)
+        self.prior = prior
+        self.mapls_chain_init = mapls_chain_init
+        self.verbose = verbose
+        self.prevalence_samples = None
+
+    def aggregation_fit(self, classif_predictions, labels):
+        self.train_post = classif_predictions
+        if self.exact_train_prev:
+            self.train_prevalence = F.prevalence_from_labels(labels, classes=self.classes_)
+        else:
+            self.train_prevalence = F.prevalence_from_probabilities(classif_predictions)
+        self.ilr = _JaxILRTransformation()
+        return self
+
+    def sample_from_posterior(self, classif_predictions):
+        n_test, n_classes = classif_predictions.shape
+        map_estimate, lam = mapls(
+            self.train_post,
+            test_probs=classif_predictions,
+            pz=self.train_prevalence,
+            return_lambda=True,
+        )
+
+        z0 = self.ilr(map_estimate)
+        if isinstance(self.prior, str) and self.prior in {'map', 'map2'}:
+            prior_context = {
+                "test_probs": classif_predictions,
+                "train_prev": self.train_prevalence,
+                "map_estimate": map_estimate,
+            }
+            alpha = _resolve_dirichlet_prior(
+                self.prior,
+                n_classes,
+                allow_mapls_priors=True,
+                n_test=n_test,
+                map_prev=prior_context,
+                map_lambda=lam,
+            )
+        else:
+            alpha = _resolve_dirichlet_prior(self.prior, n_classes)
+
+        mcmc = MCMC(
+            NUTS(self._numpyro_model),
+            num_warmup=self.num_warmup,
+            num_samples=self.num_samples,
+            num_chains=1,
+            progress_bar=self.verbose,
+        )
+        mcmc.run(
+            jrandom.PRNGKey(self.mcmc_seed),
+            test_posteriors=classif_predictions,
+            alpha=alpha,
+            init_params={"z": z0} if self.mapls_chain_init else None,
+        )
+
+        samples = mcmc.get_samples()["z"]
+        self.prevalence_samples = np.asarray(self.ilr.inverse(samples))
+        return self.prevalence_samples
+
+    def aggregate(self, classif_predictions: np.ndarray):
+        return self.sample_from_posterior(classif_predictions).mean(axis=0)
+
+    def predict_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
+        confidence_level = self.confidence_level if confidence_level is None else confidence_level
+        classif_predictions = self.classify(instances)
+        point_estimate = self.aggregate(classif_predictions)
+        region = WithConfidenceABC.construct_region(
+            self.prevalence_samples,
+            confidence_level=confidence_level,
+            method=self.region,
+        )
+        return point_estimate, region
+
+    def _log_likelihood(self, test_classif, test_prev, train_prev):
+        log_w = jnp.log(test_prev) - jnp.log(train_prev)
+        return jnp.sum(jax_logsumexp(jnp.log(test_classif) + log_w, axis=-1))
+
+    def _numpyro_model(self, test_posteriors, alpha):
+        test_posteriors = jnp.asarray(test_posteriors)
+        n_classes = test_posteriors.shape[1]
+
+        z = numpyro.sample("z", dist.Normal(jnp.zeros(n_classes - 1), 1.0))
+        prev = self.ilr.inverse(z)
+        train_prev = jnp.asarray(self.train_prevalence)
+        alpha = jnp.asarray(alpha)
+
+        numpyro.factor("dirichlet_prior", dist.Dirichlet(alpha).log_prob(prev))
+        numpyro.factor(
+            "likelihood",
+            (1.0 / self.temperature) * self._log_likelihood(test_posteriors, test_prev=prev, train_prev=train_prev),
+        )
+
+
+def mapls(
+    train_probs: np.ndarray,
+    test_probs: np.ndarray,
+    pz: np.ndarray,
+    qy_mode: str = 'soft',
+    max_iter: int = 100,
+    init_mode: str = 'identical',
+    lam: float = None,
+    dvg_name='kl',
+    return_lambda=False,
+):
+    cls_num = len(pz)
+    assert test_probs.shape[-1] == cls_num
+    if not isinstance(max_iter, int) or max_iter < 0:
+        raise ValueError(f'expected a non-negative integer for max_iter; found {max_iter!r}')
+
+    if dvg_name == 'kl':
+        dvg = kl_div
+    elif dvg_name == 'js':
+        dvg = js_div
+    else:
+        raise ValueError(f'Unsupported distribution distance measure {dvg_name!r}; expected "kl" or "js"')
+
+    q_prior = np.ones(cls_num) / cls_num
+    if lam is None:
+        lam = get_lambda(test_probs, pz, q_prior, dvg=dvg, max_iter=max_iter)
+
+    qz = mapls_em(
+        test_probs,
+        pz,
+        lam,
+        q_prior,
+        cls_num,
+        init_mode=init_mode,
+        max_iter=max_iter,
+        qy_mode=qy_mode,
+    )
+    return (qz, lam) if return_lambda else qz
+
+
+def mapls_em(probs, pz, lam, q_prior, cls_num, init_mode='identical', max_iter=100, qy_mode='soft'):
+    pz = np.asarray(pz, dtype=float)
+    pz = pz / np.sum(pz)
+    if init_mode == 'uniform':
+        qz = np.ones(cls_num) / cls_num
+    elif init_mode == 'identical':
+        qz = pz.copy()
+    else:
+        raise ValueError('init_mode should be either "uniform" or "identical"')
+
+    w = qz / pz
+    for _ in range(max_iter):
+        mapls_probs = normalized(probs * w, axis=-1, order=1)
+        if qy_mode == 'hard':
+            pred = np.argmax(mapls_probs, axis=-1)
+            qz_new = np.bincount(pred.reshape(-1), minlength=cls_num)
+        elif qy_mode == 'soft':
+            qz_new = np.mean(mapls_probs, axis=0)
+        else:
+            raise ValueError('qy_mode should be either "soft" or "hard"')
+
+        qz = lam * qz_new + (1 - lam) * q_prior
+        qz /= qz.sum()
+        w = qz / pz
+
+    return qz
+
+
+def get_lambda(test_probs, pz, q_prior, dvg, max_iter=50):
+    n_classes = len(pz)
+    qz_pred = mapls_em(test_probs, pz, 1, 0, n_classes, max_iter=max_iter)
+
+    tu_div = dvg(qz_pred, q_prior)
+    ts_div = dvg(qz_pred, pz)
+    su_div = dvg(pz, q_prior)
+
+    su_conf = 1 - lambda_forward(su_div, lambda_inverse(dpq=0.5, lam=0.2))
+    tu_conf = lambda_forward(tu_div, lambda_inverse(dpq=0.5, lam=su_conf))
+    ts_conf = lambda_forward(ts_div, lambda_inverse(dpq=0.5, lam=su_conf))
+
+    confs = np.array([tu_conf, 1 - ts_conf])
+    weights = np.array([0.9, 0.1])
+    return np.sum(weights * confs)
+
+
+def lambda_inverse(dpq, lam):
+    return (1 / (1 - lam) - 1) / dpq
+
+
+def lambda_forward(dpq, gamma):
+    return gamma * dpq / (1 + gamma * dpq)
+
+
+def get_lamda(test_probs, pz, q_prior, dvg, max_iter=50):
+    return get_lambda(test_probs, pz, q_prior, dvg, max_iter=max_iter)
+
+
+def lam_inv(dpq, lam):
+    return lambda_inverse(dpq, lam)
+
+
+def lam_forward(dpq, gamma):
+    return lambda_forward(dpq, gamma)
+
+
+def kl_div(p, q, eps=1e-12):
+    p = np.asarray(p, dtype=float)
+    q = np.asarray(q, dtype=float)
+
+    mask = p > 0
+    return np.sum(p[mask] * np.log(p[mask] / (q[mask] + eps)))
+
+
+def js_div(p, q):
+    assert (np.abs(np.sum(p) - 1) < 1e-6) and (np.abs(np.sum(q) - 1) < 1e-6)
+    m = (p + q) / 2
+    return kl_div(p, m) / 2 + kl_div(q, m) / 2
+
+
+def normalized(a, axis=-1, order=2):
+    l2 = np.atleast_1d(np.linalg.norm(a, order, axis))
+    l2[l2 == 0] = 1
+    return a / np.expand_dims(l2, axis)
+
+
+def alpha0_from_lambda(lam, n_test, n_classes):
+    return 1 + n_test * (1 - lam) / (lam * n_classes)
+
+
+def alpha0_from_lamda(lam, n_test, n_classes):
+    return alpha0_from_lambda(lam, n_test, n_classes)
+
+
+def calibrate_temperature(
+    method: WithConfidenceABC,
+    train: LabelledCollection,
+    val_prot: AbstractProtocol,
+    temp_grid=(0.5, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0),
+    nominal_coverage: float = 0.95,
+    amplitude_threshold=1.0,
+    criterion: str = 'winkler',
+    n_jobs: int = 1,
+    verbose: bool = True,
+):
+    """
+    Calibrates the temperature parameter of a Bayesian quantifier with
+    confidence regions by selecting the value that yields the best validation
+    trade-off between nominal coverage and region sharpness.
+
+    The method is first fitted on ``train``. For each candidate temperature,
+    the fitted quantifier is deep-copied, its ``temperature`` attribute is
+    replaced, and it is evaluated on the samples generated by ``val_prot``.
+    Candidate temperatures whose average region amplitude exceeds
+    ``amplitude_threshold`` are discarded.
+
+    When ``criterion='winkler'``, the surviving candidate with minimum mean
+    Winkler score is selected. When ``criterion='auto'``, the selected
+    temperature is the one whose empirical coverage is closest to
+    ``nominal_coverage``.
+
+    :param method: a quantifier implementing :class:`WithConfidenceABC` and
+        exposing a writable ``temperature`` attribute
+    :param train: training set used to fit the quantifier
+    :param val_prot: validation protocol yielding pairs ``(sample, true_prev)``
+    :param temp_grid: candidate temperatures to evaluate
+    :param nominal_coverage: target confidence level used by the Winkler score
+        and coverage selection
+    :param amplitude_threshold: maximum allowed average simplex proportion of
+        the region. It can also be set to ``'auto'`` to use a heuristic based
+        on the number of classes
+    :param criterion: either ``'winkler'`` (default) or ``'auto'``
+    :param n_jobs: number of parallel jobs across candidate temperatures
+    :param verbose: whether to display progress information
+    :return: the selected temperature value
+    """
+    if not hasattr(method, 'temperature'):
+        raise ValueError(f'{method.__class__.__name__} does not expose a temperature attribute')
+    if not isinstance(method, WithConfidenceABC):
+        raise TypeError(f'{method.__class__.__name__} is not an instance of WithConfidenceABC')
+    if not 0 < nominal_coverage < 1:
+        raise ValueError(f'{nominal_coverage=} must be in the interval (0,1)')
+    if criterion not in {'auto', 'winkler'}:
+        raise ValueError(f'unknown {criterion=}; valid ones are "auto" or "winkler"')
+    if amplitude_threshold != 'auto':
+        if not isinstance(amplitude_threshold, Real) or amplitude_threshold > 1.0:
+            raise ValueError(
+                f'wrong value for {amplitude_threshold=}; it must either be "auto" or a real value <= 1.0'
+            )
+    temperatures = sorted(_validate_temperature(temp) for temp in temp_grid)
+
+    if amplitude_threshold == 'auto':
+        amplitude_threshold = 0.1 / np.log(train.n_classes + 1)
+
+    if amplitude_threshold > 0.1:
+        print(f'warning: the {amplitude_threshold=} is too large; this may lead to uninformative regions')
+
+    def _evaluate_temperature_job(job_id, temperature):
+        local_method = copy.deepcopy(method)
+        local_method.temperature = temperature
+
+        coverage = 0
+        amplitudes = []
+        winklers = []
+        errors = []
+
+        pbar = tqdm(
+            enumerate(val_prot()),
+            position=job_id,
+            total=val_prot.total(),
+            disable=not verbose,
+        )
+
+        for i, (sample, prev) in pbar:
+            point_estim, conf_region = local_method.predict_conf(sample)
+
+            if prev in conf_region:
+                coverage += 1
+
+            amplitudes.append(conf_region.montecarlo_proportion(n_trials=50_000))
+            if criterion == 'winkler':
+                winklers.append(conf_region.mean_winkler_score(true_prev=prev, alpha=1 - nominal_coverage))
+            errors.append(qp.error.mae(prev, point_estim))
+
+            description = (
+                f'job={job_id} T={temperature}: '
+                f'MAE={np.mean(errors):.6f} '
+                f'coverage={coverage / (i + 1) * 100:.2f}% '
+                f'amplitude={np.mean(amplitudes) * 100:.4f}% '
+            )
+            if criterion == 'winkler':
+                description += f'winkler={np.mean(winklers):.4f}'
+            pbar.set_description(description)
+
+        mean_coverage = coverage / val_prot.total()
+        mean_amplitude = float(np.mean(amplitudes))
+        mean_winkler = float(np.mean(winklers)) if criterion == 'winkler' else None
+        mean_error = float(np.mean(errors))
+        return temperature, mean_coverage, mean_amplitude, mean_winkler, mean_error
+
+    method.fit(*train.Xy)
+    raw_results = Parallel(n_jobs=n_jobs, backend="loky")(
+        delayed(_evaluate_temperature_job)(job_id, temperature)
+        for job_id, temperature in tqdm(enumerate(temperatures), disable=not verbose)
+    )
+    filtered_results = [
+        (temperature, coverage, amplitude, winkler, error)
+        for temperature, coverage, amplitude, winkler, error in raw_results
+        if amplitude < amplitude_threshold
+    ]
+
+    chosen_temperature = 1.0
+    chosen_coverage = chosen_amplitude = chosen_winkler = chosen_error = None
+
+    if filtered_results:
+        if criterion == 'winkler':
+            chosen_temperature, chosen_coverage, chosen_amplitude, chosen_winkler, chosen_error = min(
+                filtered_results, key=lambda item: item[3]
+            )
+        else:
+            chosen_temperature, chosen_coverage, chosen_amplitude, chosen_winkler, chosen_error = min(
+                filtered_results, key=lambda item: abs(item[1] - nominal_coverage)
+            )
+
+    if verbose and chosen_coverage is not None:
+        message = (
+            f'\nChosen_temperature={chosen_temperature:.2f} got '
+            f'MAE={chosen_error:.6f} '
+            f'coverage={chosen_coverage * 100:.2f}% '
+            f'amplitude={chosen_amplitude * 100:.4f}% '
+        )
+        if criterion == 'winkler':
+            message += f'winkler={chosen_winkler:.4f}'
+        print(message)
+
+    return chosen_temperature
+
+
+def temp_calibration(*args, **kwargs):
+    """
+    Backward-compatible alias for :func:`calibrate_temperature`.
+    """
+    return calibrate_temperature(*args, **kwargs)

quapy/method/_kdey.py

@@ -1,11 +1,12 @@
 import numpy as np
+from numbers import Real
 from sklearn.base import BaseEstimator
 from sklearn.neighbors import KernelDensity
 
 import quapy as qp
 from quapy.method.aggregative import AggregativeSoftQuantifier
 import quapy.functional as F
-
+from scipy.special import logsumexp
 from sklearn.metrics.pairwise import rbf_kernel
 
 
@@ -15,44 +16,57 @@ class KDEBase:
     """
 
     BANDWIDTH_METHOD = ['scott', 'silverman']
+    KERNELS = ['gaussian', 'aitchison', 'ilr']
 
     @classmethod
-    def _check_bandwidth(cls, bandwidth):
+    def _check_bandwidth(cls, bandwidth, kernel):
         """
         Checks that the bandwidth parameter is correct
 
         :param bandwidth: either a string (see BANDWIDTH_METHOD) or a float
         :return: the bandwidth if the check is passed, or raises an exception for invalid values
         """
-        assert bandwidth in KDEBase.BANDWIDTH_METHOD or isinstance(bandwidth, float), \
+        assert bandwidth in KDEBase.BANDWIDTH_METHOD or isinstance(bandwidth, Real), \
             f'invalid bandwidth, valid ones are {KDEBase.BANDWIDTH_METHOD} or float values'
-        if isinstance(bandwidth, float):
-            assert 0 < bandwidth < 1,  \
-                "the bandwith for KDEy should be in (0,1), since this method models the unit simplex"
+        if isinstance(bandwidth, Real):
+            bandwidth = float(bandwidth)
         return bandwidth
 
-    def get_kde_function(self, X, bandwidth):
+    @classmethod
+    def _check_kernel(cls, kernel):
+        assert kernel in KDEBase.KERNELS, f'unknown {kernel=}'
+        return kernel
+
+    def get_kde_function(self, X, bandwidth, kernel):
         """
         Wraps the KDE function from scikit-learn.
 
         :param X: data for which the density function is to be estimated
         :param bandwidth: the bandwidth of the kernel
+        :param kernel: the kernel family
         :return: a scikit-learn's KernelDensity object
         """
+        X = self.transform_posteriors(X, kernel)
+        bandwidth = self.effective_bandwidth(bandwidth, kernel)
         return KernelDensity(bandwidth=bandwidth).fit(X)
 
-    def pdf(self, kde, X):
+    def pdf(self, kde, X, kernel, log_densities=False):
         """
         Wraps the density evalution of scikit-learn's KDE. Scikit-learn returns log-scores (s), so this
         function returns :math:`e^{s}`
 
         :param kde: a previously fit KDE function
         :param X: the data for which the density is to be estimated
+        :param kernel: the kernel family
         :return: np.ndarray with the densities
         """
-        return np.exp(kde.score_samples(X))
+        X = self.transform_posteriors(X, kernel)
+        log_density = kde.score_samples(X)
+        if log_densities:
+            return log_density
+        return np.exp(log_density)
 
-    def get_mixture_components(self, X, y, classes, bandwidth):
+    def get_mixture_components(self, X, y, classes, bandwidth, kernel):
         """
         Returns an array containing the mixture components, i.e., the KDE functions for each class.
 
@@ -60,15 +74,50 @@ class KDEBase:
         :param y: the class labels
         :param n_classes: integer, the number of classes
         :param bandwidth: float, the bandwidth of the kernel
+        :param kernel: the kernel family
         :return: a list of KernelDensity objects, each fitted with the corresponding class-specific covariates
         """
         class_cond_X = []
         for cat in classes:
             selX = X[y==cat]
             if selX.size==0:
-                selX = [F.uniform_prevalence(len(classes))]
+                raise ValueError(f'empty class {cat}')
             class_cond_X.append(selX)
-        return [self.get_kde_function(X_cond_yi, bandwidth) for X_cond_yi in class_cond_X]
+        return [self.get_kde_function(X_cond_yi, bandwidth, kernel) for X_cond_yi in class_cond_X]
+
+    def transform_posteriors(self, X, kernel):
+        if kernel in {'aitchison', 'ilr'}:
+            X = self.shrink_posteriors(X)
+        if kernel == 'aitchison':
+            return self.clr_transform(X)
+        if kernel == 'ilr':
+            return self.ilr_transform(X)
+        return X
+
+    def shrink_posteriors(self, X):
+        shrinkage = getattr(self, 'shrinkage', 0.0)
+        if shrinkage <= 0:
+            return X
+        X = np.asarray(X)
+        n_classes = X.shape[-1]
+        uniform = np.full(n_classes, 1.0 / n_classes, dtype=X.dtype)
+        return (1.0 - shrinkage) * X + shrinkage * uniform
+
+    def effective_bandwidth(self, bandwidth, kernel):
+        shrinkage = getattr(self, 'shrinkage', 0.0)
+        if shrinkage > 0 and kernel in {'aitchison', 'ilr'} and isinstance(bandwidth, Real):
+            return (1.0 - shrinkage) * float(bandwidth)
+        return bandwidth
+
+    def clr_transform(self, X):
+        if not hasattr(self, 'clr'):
+            self.clr = F.CLRtransformation()
+        return self.clr(X)
+
+    def ilr_transform(self, X):
+        if not hasattr(self, 'ilr'):
+            self.ilr = F.ILRtransformation()
+        return self.ilr(X)
 
 
 class KDEyML(AggregativeSoftQuantifier, KDEBase):
@@ -107,17 +156,31 @@ class KDEyML(AggregativeSoftQuantifier, KDEBase):
         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 bandwidth: float, the bandwidth of the Kernel
+    :param kernel: kernel of KDE, valid ones are in KDEBase.KERNELS
+    :param shrinkage: amount of shrinkage towards the uniform distribution to apply before
+        Aitchison/ILR transformations. Must be in ``[0,1)``.
     :param random_state: a seed to be set before fitting any base quantifier (default None)
     """
 
     def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5, bandwidth=0.1,
-                 random_state=None):
+                 kernel='gaussian', shrinkage=0.0, random_state=None):
         super().__init__(classifier, fit_classifier, val_split)
-        self.bandwidth = KDEBase._check_bandwidth(bandwidth)
+        self.bandwidth = KDEBase._check_bandwidth(bandwidth, kernel)
+        self.kernel = self._check_kernel(kernel)
+        assert 0 <= shrinkage < 1, 'shrinkage must be in [0,1)'
+        assert self.kernel != 'gaussian' or shrinkage == 0, \
+            'shrinkage is only supported for Aitchison/ILR kernels'
+        self.shrinkage = float(shrinkage)
         self.random_state=random_state
 
     def aggregation_fit(self, classif_predictions, labels):
-        self.mix_densities = self.get_mixture_components(classif_predictions, labels, self.classes_, self.bandwidth)
+        self.mix_densities = self.get_mixture_components(
+            classif_predictions,
+            labels,
+            self.classes_,
+            self.bandwidth,
+            self.kernel,
+        )
         return self
 
     def aggregate(self, posteriors: np.ndarray):
@@ -129,15 +192,25 @@ class KDEyML(AggregativeSoftQuantifier, KDEBase):
         :return: a vector of class prevalence estimates
         """
         with qp.util.temp_seed(self.random_state):
-            epsilon = 1e-10
+            epsilon = 1e-12
             n_classes = len(self.mix_densities)
-            test_densities = [self.pdf(kde_i, posteriors) for kde_i in self.mix_densities]
+            if (self.kernel != 'gaussian' and n_classes >= 20) or n_classes >= 30:
+                test_log_densities = [
+                    self.pdf(kde_i, posteriors, self.kernel, log_densities=True)
+                    for kde_i in self.mix_densities
+                ]
 
-            def neg_loglikelihood(prev):
-                # test_mixture_likelihood = sum(prev_i * dens_i for prev_i, dens_i in zip (prev, test_densities))
-                test_mixture_likelihood = prev @ test_densities
-                test_loglikelihood = np.log(test_mixture_likelihood + epsilon)
-                return  -np.sum(test_loglikelihood)
+                def neg_loglikelihood(prev):
+                    prev = qp.error.smooth(prev, eps=epsilon)
+                    test_loglikelihood = logsumexp(np.log(prev)[:, None] + test_log_densities, axis=0)
+                    return -np.sum(test_loglikelihood)
+            else:
+                test_densities = [self.pdf(kde_i, posteriors, self.kernel) for kde_i in self.mix_densities]
+
+                def neg_loglikelihood(prev):
+                    test_mixture_likelihood = prev @ test_densities
+                    test_loglikelihood = np.log(test_mixture_likelihood + epsilon)
+                    return -np.sum(test_loglikelihood)
 
             return F.optim_minimize(neg_loglikelihood, n_classes)
 
@@ -192,18 +265,22 @@ class KDEyHD(AggregativeSoftQuantifier, KDEBase):
 
         super().__init__(classifier, fit_classifier, val_split)
         self.divergence = divergence
-        self.bandwidth = KDEBase._check_bandwidth(bandwidth)
+        self.bandwidth = KDEBase._check_bandwidth(bandwidth, kernel='gaussian')
         self.random_state=random_state
         self.montecarlo_trials = montecarlo_trials
 
     def aggregation_fit(self, classif_predictions, labels):
-        self.mix_densities = self.get_mixture_components(classif_predictions, labels, self.classes_, self.bandwidth)
+        self.mix_densities = self.get_mixture_components(
+            classif_predictions, labels, self.classes_, self.bandwidth, 'gaussian'
+        )
 
         N = self.montecarlo_trials
         rs = self.random_state
         n = len(self.classes_)
         self.reference_samples = np.vstack([kde_i.sample(N//n, random_state=rs) for kde_i in self.mix_densities])
-        self.reference_classwise_densities = np.asarray([self.pdf(kde_j, self.reference_samples) for kde_j in self.mix_densities])
+        self.reference_classwise_densities = np.asarray(
+            [self.pdf(kde_j, self.reference_samples, 'gaussian') for kde_j in self.mix_densities]
+        )
         self.reference_density = np.mean(self.reference_classwise_densities, axis=0)  # equiv. to (uniform @ self.reference_classwise_densities)
 
         return self
@@ -213,8 +290,8 @@ class KDEyHD(AggregativeSoftQuantifier, KDEBase):
         # apply importance sampling (IS). In this version we compute D(p_alpha||q) with IS
         n_classes = len(self.mix_densities)
 
-        test_kde = self.get_kde_function(posteriors, self.bandwidth)
-        test_densities = self.pdf(test_kde, self.reference_samples)
+        test_kde = self.get_kde_function(posteriors, self.bandwidth, 'gaussian')
+        test_densities = self.pdf(test_kde, self.reference_samples, 'gaussian')
 
         def f_squared_hellinger(u):
             return (np.sqrt(u)-1)**2
@@ -279,7 +356,7 @@ class KDEyCS(AggregativeSoftQuantifier):
 
     def __init__(self, classifier: BaseEstimator=None, fit_classifier=True, val_split=5, bandwidth=0.1):
         super().__init__(classifier, fit_classifier, val_split)
-        self.bandwidth = KDEBase._check_bandwidth(bandwidth)
+        self.bandwidth = KDEBase._check_bandwidth(bandwidth, kernel='gaussian')
 
     def gram_matrix_mix_sum(self, X, Y=None):
         # this adapts the output of the rbf_kernel function (pairwise evaluations of Gaussian kernels k(x,y))
@@ -354,4 +431,3 @@ class KDEyCS(AggregativeSoftQuantifier):
             return partA + partB #+ partC
 
         return F.optim_minimize(divergence, n)
-

quapy/method/aggregative.py

@@ -3,7 +3,6 @@ from argparse import ArgumentError
 from copy import deepcopy
 from typing import Callable, Literal, Union
 import numpy as np
-from abstention.calibration import NoBiasVectorScaling, TempScaling, VectorScaling
 from numpy.f2py.crackfortran import true_intent_list
 from sklearn.base import BaseEstimator
 from sklearn.calibration import CalibratedClassifierCV
@@ -18,13 +17,27 @@ from quapy.functional import get_divergence
 from quapy.classification.svmperf import SVMperf
 from quapy.data import LabelledCollection
 from quapy.method.base import BaseQuantifier, BinaryQuantifier, OneVsAllGeneric
-from quapy.method import _bayesian
 
 # import warnings
 # from sklearn.exceptions import ConvergenceWarning
 # warnings.filterwarnings("ignore", category=ConvergenceWarning)
 
 
+def _get_abstention_calibrators():
+    try:
+        from abstention.calibration import NoBiasVectorScaling, TempScaling, VectorScaling
+    except ImportError as exc:
+        raise ImportError(
+            "Posterior calibration for EMQ requires the optional 'abstention' package."
+        ) from exc
+    return {
+        'nbvs': NoBiasVectorScaling(),
+        'bcts': TempScaling(bias_positions='all'),
+        'ts': TempScaling(),
+        'vs': VectorScaling(),
+    }
+
+
 # Abstract classes
 # ------------------------------------
 
@@ -849,12 +862,7 @@ class EMQ(AggregativeSoftQuantifier):
                                 "validation data")
 
         if self.calib is not None:
-            calibrator = {
-                'nbvs': NoBiasVectorScaling(),
-                'bcts': TempScaling(bias_positions='all'),
-                'ts': TempScaling(),
-                'vs': VectorScaling()
-            }.get(self.calib, None)
+            calibrator = _get_abstention_calibrators().get(self.calib, None)
 
             if calibrator is None:
                 raise ValueError(f'invalid value for {self.calib=}; valid ones are {EMQ.CALIB_OPTIONS}')

quapy/method/confidence.py

@@ -1,24 +1,34 @@
+from numbers import Number
+from typing import Iterable
+
 import numpy as np
+from joblib import Parallel, delayed
 from sklearn.base import BaseEstimator
 from sklearn.metrics import confusion_matrix
 
 import quapy as qp
 import quapy.functional as F
-from quapy.method import _bayesian
+from quapy.functional import CompositionalTransformation, CLRtransformation, ILRtransformation
 from quapy.data import LabelledCollection
 from quapy.method.aggregative import AggregativeQuantifier, AggregativeCrispQuantifier, AggregativeSoftQuantifier, BinaryAggregativeQuantifier
 from scipy.stats import chi2
 from sklearn.utils import resample
 from abc import ABC, abstractmethod
-from scipy.special import softmax, factorial
+from scipy.special import factorial
 import copy
 from functools import lru_cache
+from tqdm import tqdm
 
 """
 This module provides implementation of different types of confidence regions, and the implementation of Bootstrap
 for AggregativeQuantifiers.
 """
 
+
+def _get_bayesian_module():
+    from quapy.method import _bayesian
+    return _bayesian
+
 class ConfidenceRegionABC(ABC):
     """
     Abstract class of confidence regions
@@ -85,12 +95,58 @@ class ConfidenceRegionABC(ABC):
         """ Returns internal samples """
         ...
 
+    def __contains__(self, p):
+        """
+        Overloads in operator, checks if `p` is contained in the region
+
+        :param p: array-like
+        :return: boolean
+        """
+        p = np.asarray(p)
+        assert p.ndim==1, f'unexpected shape for point parameter'
+        return self.coverage(p)==1.
+
+    def closest_point_in_region(self, p, tol=1e-6, max_iter=30):
+        """
+        Finds the closes point to p that belongs to the region. Assumes the region is convex.
+
+        :param p: array-like, the point
+        :param tol: float, error tolerance
+        :param max_iter: int, max number of iterations
+        :returns: array-like, the closes point to p in the segment between p and the center of the region, that
+            belongs to the region
+        """
+        p = np.asarray(p, dtype=float)
+
+        # if p in region, returns p itself
+        if p in self:
+            return p.copy()
+
+        # center of the region
+        c = self.point_estimate()
+
+        # binary search in [0,1], interpolation parameter
+        # low=closest to p, high=closest to c
+        low, high = 0.0, 1.0
+        for _ in range(max_iter):
+            mid = 0.5 * (low + high)
+            x = p*(1-mid) + c*mid
+            if x in self:
+                high = mid
+            else:
+                low = mid
+            if high - low < tol:
+                break
+
+        in_boundary = p*(1-high) + c*high
+        return in_boundary
+
 
 class WithConfidenceABC(ABC):
     """
     Abstract class for confidence regions.
     """
-    METHODS = ['intervals', 'ellipse', 'ellipse-clr']
+    REGION_TYPE = ['intervals', 'ellipse', 'ellipse-clr', 'ellipse-ilr']
 
     @abstractmethod
     def predict_conf(self, instances, confidence_level=0.95) -> (np.ndarray, ConfidenceRegionABC):
@@ -118,7 +174,7 @@ class WithConfidenceABC(ABC):
         return self.predict_conf(instances=instances, confidence_level=confidence_level)
 
     @classmethod
-    def construct_region(cls, prev_estims, confidence_level=0.95, method='intervals'):
+    def construct_region(cls, prev_estims, confidence_level=0.95, method='intervals')->ConfidenceRegionABC:
         """
         Construct a confidence region given many prevalence estimations.
 
@@ -136,6 +192,8 @@ class WithConfidenceABC(ABC):
             region = ConfidenceEllipseSimplex(prev_estims, confidence_level=confidence_level)
         elif method == 'ellipse-clr':
             region = ConfidenceEllipseCLR(prev_estims, confidence_level=confidence_level)
+        elif method == 'ellipse-ilr':
+            region = ConfidenceEllipseILR(prev_estims, confidence_level=confidence_level)
 
         if region is None:
             raise NotImplementedError(f'unknown method {method}')
@@ -153,7 +211,7 @@ def simplex_volume(n):
     return 1 / factorial(n)
 
 
-def within_ellipse_prop(values, mean, prec_matrix, chi2_critical):
+def within_ellipse_prop__(values, mean, prec_matrix, chi2_critical):
     """
     Checks the proportion of values that belong to the ellipse with center `mean` and precision matrix `prec_matrix`
     at a distance `chi2_critical`.
@@ -186,102 +244,91 @@ def within_ellipse_prop(values, mean, prec_matrix, chi2_critical):
     return within_elipse * 1.0
 
 
-class ConfidenceEllipseSimplex(ConfidenceRegionABC):
+def within_ellipse_prop(values, mean, prec_matrix, chi2_critical):
     """
-    Instantiates a Confidence Ellipse in the probability simplex.
+        Checks the proportion of values that belong to the ellipse with center `mean` and precision matrix `prec_matrix`
+        at a distance `chi2_critical`.
 
-    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
-    :param confidence_level: float, the confidence level (default 0.95)
+        :param values: a np.ndarray of shape (n_dim,) or (n_values, n_dim,)
+        :param mean: a np.ndarray of shape (n_dim,) with the center of the ellipse
+        :param prec_matrix: a np.ndarray with the precision matrix (inverse of the
+            covariance matrix) of the ellipse. If this inverse cannot be computed
+            then None must be passed
+        :param chi2_critical: float, the chi2 critical value
+
+        :return: float in [0,1], the fraction of values that are contained in the ellipse
+            defined by the mean (center), the precision matrix (shape), and the chi2_critical value (distance).
+            If `values` is only one value, then either 0. (not contained) or 1. (contained) is returned.
+        """
+    if prec_matrix is None:
+        return 0.
+
+    values = np.atleast_2d(values)
+    diff = values - mean
+    d_M_squared = np.sum(diff @ prec_matrix * diff, axis=-1)
+    within_ellipse = d_M_squared <= chi2_critical
+
+    if len(within_ellipse) == 1:
+        return float(within_ellipse[0])
+    else:
+        return float(np.mean(within_ellipse))
+
+
+def closest_point_on_ellipsoid(p, mean, cov, chi2_critical, tol=1e-9, max_iter=100):
+    """
+    Computes the closest point on the ellipsoid defined by:
+        (x - mean)^T cov^{-1} (x - mean) = chi2_critical
     """
 
-    def __init__(self, samples, confidence_level=0.95):
+    p = np.asarray(p)
+    mean = np.asarray(mean)
+    Sigma = np.asarray(cov)
 
-        assert 0. < confidence_level < 1., f'{confidence_level=} must be in range(0,1)'
+    # Precompute precision matrix
+    P = np.linalg.pinv(Sigma)
+    d = P.shape[0]
 
-        samples = np.asarray(samples)
+    # Define v = p - mean
+    v = p - mean
 
-        self.mean_ = samples.mean(axis=0)
-        self.cov_ = np.cov(samples, rowvar=False, ddof=1)
+    # If p is inside the ellipsoid, return p itself
+    M_dist = v @ P @ v
+    if M_dist <= chi2_critical:
+        return p.copy()
 
-        try:
-            self.precision_matrix_ = np.linalg.inv(self.cov_)
-        except:
-            self.precision_matrix_ = None
+    # Function to compute x(lambda)
+    def x_lambda(lmbda):
+        A = np.eye(d) + lmbda * P
+        return mean + np.linalg.solve(A, v)
 
-        self.dim = samples.shape[-1]
-        self.ddof = self.dim - 1
+    # Function whose root we want: f(lambda) = Mahalanobis distance - chi2
+    def f(lmbda):
+        x = x_lambda(lmbda)
+        diff = x - mean
+        return diff @ P @ diff - chi2_critical
 
-        # critical chi-square value
-        self.confidence_level = confidence_level
-        self.chi2_critical_ = chi2.ppf(confidence_level, df=self.ddof)
-        self._samples = samples
+    # Bisection search over lambda >= 0
+    l_low, l_high = 0.0, 1.0
 
-    @property
-    def samples(self):
-        return self._samples
+    # Increase high until f(high) < 0
+    while f(l_high) > 0:
+        l_high *= 2
+        if l_high > 1e12:
+            raise RuntimeError("Failed to bracket the root.")
 
-    def point_estimate(self):
-        """
-        Returns the point estimate, the center of the ellipse.
+    # Bisection
+    for _ in range(max_iter):
+        l_mid = 0.5 * (l_low + l_high)
+        fm = f(l_mid)
+        if abs(fm) < tol:
+            break
+        if fm > 0:
+            l_low = l_mid
+        else:
+            l_high = l_mid
 
-        :return: np.ndarray of shape (n_classes,)
-        """
-        return self.mean_
-
-    def coverage(self, true_value):
-        """
-        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
-        fraction of these that are contained in the region, if more than one value is passed. If only one value is
-        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
-
-        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
-        :return: float in [0,1]
-        """
-        return within_ellipse_prop(true_value, self.mean_, self.precision_matrix_, self.chi2_critical_)
-
-
-class ConfidenceEllipseCLR(ConfidenceRegionABC):
-    """
-    Instantiates a Confidence Ellipse in the Centered-Log Ratio (CLR) space.
-
-    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
-    :param confidence_level: float, the confidence level (default 0.95)
-    """
-
-    def __init__(self, samples, confidence_level=0.95):
-        samples = np.asarray(samples)
-        self.clr = CLRtransformation()
-        Z = self.clr(samples)
-        self.mean_ = np.mean(samples, axis=0)
-        self.conf_region_clr = ConfidenceEllipseSimplex(Z, confidence_level=confidence_level)
-        self._samples = samples
-
-    @property
-    def samples(self):
-        return self._samples
-
-    def point_estimate(self):
-        """
-        Returns the point estimate, the center of the ellipse.
-
-        :return: np.ndarray of shape (n_classes,)
-        """
-        # The inverse of the CLR does not coincide with the true mean, because the geometric mean
-        # requires smoothing the prevalence vectors and this affects the softmax (inverse);
-        # return self.clr.inverse(self.mean_) # <- does not coincide
-        return self.mean_
-
-    def coverage(self, true_value):
-        """
-        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
-        fraction of these that are contained in the region, if more than one value is passed. If only one value is
-        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
-
-        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
-        :return: float in [0,1]
-        """
-        transformed_values = self.clr(true_value)
-        return self.conf_region_clr.coverage(transformed_values)
+    l_opt = l_mid
+    return x_lambda(l_opt)
 
 
 class ConfidenceIntervals(ConfidenceRegionABC):
@@ -290,18 +337,30 @@ class ConfidenceIntervals(ConfidenceRegionABC):
 
     :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
     :param confidence_level: float, the confidence level (default 0.95)
+    :param bonferroni_correction: bool (default False), if True, a Bonferroni correction
+        is applied to the significance level (`alpha`) before computing confidence intervals.
+        The correction consists of replacing `alpha` with `alpha/n_classes`. When
+        `n_classes=2` the correction is not applied because there is only one verification test
+        since the other class is constrained. This is not necessarily true for n_classes>2.
     """
-    def __init__(self, samples, confidence_level=0.95):
+    def __init__(self, samples, confidence_level=0.95, bonferroni_correction=False):
         assert 0 < confidence_level < 1, f'{confidence_level=} must be in range(0,1)'
+        assert samples.ndim == 2, 'unexpected shape; must be (n_bootstrap_samples, n_classes)'
 
         samples = np.asarray(samples)
 
         self.means_ = samples.mean(axis=0)
+        self.confidence_level = confidence_level
         alpha = 1-confidence_level
+        if bonferroni_correction:
+            n_classes = samples.shape[-1]
+            if n_classes>2:
+                alpha = alpha/n_classes
         low_perc = (alpha/2.)*100
         high_perc = (1-alpha/2.)*100
         self.I_low, self.I_high = np.percentile(samples, q=[low_perc, high_perc], axis=0)
         self._samples = samples
+        self.alpha = alpha
 
     @property
     def samples(self):
@@ -330,36 +389,278 @@ class ConfidenceIntervals(ConfidenceRegionABC):
 
         return proportion
 
+    def coverage_soft(self, true_value):
+        within_intervals = np.logical_and(self.I_low <= true_value, true_value <= self.I_high)
+        return np.mean(within_intervals.astype(float))
+
     def __repr__(self):
         return '['+', '.join(f'({low:.4f}, {high:.4f})' for (low,high) in zip(self.I_low, self.I_high))+']'
 
+    @property
+    def n_dim(self):
+        return len(self.I_low)
 
-class CLRtransformation:
+    def winkler_scores(self, true_prev, alpha=None, add_ae=False):
+        true_prev = np.asarray(true_prev)
+        assert true_prev.ndim == 1, 'unexpected dimensionality for true_prev'
+        assert len(true_prev)==self.n_dim, \
+            f'unexpected number of dimensions; found {true_prev.ndim}, expected {self.n_dim}'
+
+        def winkler_score(low, high, true_val, alpha, center):
+            amp = high-low
+            scale_cost = 2./alpha
+            cost = np.max([0, low-true_val], axis=0) + np.max([0, true_val-high], axis=0)
+            err = 0
+            if add_ae:
+                err = abs(true_val - center)
+            return amp + scale_cost*cost + err
+
+        alpha = alpha or self.alpha
+        return np.asarray(
+            [winkler_score(low_i, high_i, true_v, alpha, center)
+                for (low_i, high_i, true_v, center) in zip(self.I_low, self.I_high, true_prev, self.point_estimate())]
+        )
+
+    def mean_winkler_score(self, true_prev, alpha=None, add_ae=False):
+        return np.mean(self.winkler_scores(true_prev, alpha=alpha, add_ae=add_ae))
+
+
+
+class ConfidenceEllipseSimplex(ConfidenceRegionABC):
     """
-    Centered log-ratio, from component analysis
+    Instantiates a Confidence Ellipse in the probability simplex.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
     """
-    def __call__(self, X, epsilon=1e-6):
-        """
-        Applies the CLR function to X thus mapping the instances, which are contained in `\\mathcal{R}^{n}` but
-        actually lie on a `\\mathcal{R}^{n-1}` simplex, onto an unrestricted space in :math:`\\mathcal{R}^{n}`
 
-        :param X: np.ndarray of (n_instances, n_dimensions) to be transformed
-        :param epsilon: small float for prevalence smoothing
-        :return: np.ndarray of (n_instances, n_dimensions), the CLR-transformed points
-        """
-        X = np.asarray(X)
-        X = qp.error.smooth(X, epsilon)
-        G = np.exp(np.mean(np.log(X), axis=-1, keepdims=True))  # geometric mean
-        return np.log(X / G)
+    def __init__(self, samples, confidence_level=0.95):
 
-    def inverse(self, X):
-        """
-        Inverse function. However, clr.inverse(clr(X)) does not exactly coincide with X due to smoothing.
+        assert 0. < confidence_level < 1., f'{confidence_level=} must be in range(0,1)'
 
-        :param X: np.ndarray of (n_instances, n_dimensions) to be transformed
-        :return: np.ndarray of (n_instances, n_dimensions), the CLR-transformed points
+        samples = np.asarray(samples)
+
+        self.confidence_level = confidence_level
+        self.mean_ = samples.mean(axis=0)
+        self.cov_ = np.cov(samples, rowvar=False, ddof=1)
+
+        try:
+            self.precision_matrix_ = np.linalg.pinv(self.cov_)
+        except:
+            self.precision_matrix_ = None
+
+        self.dim = samples.shape[-1]
+        self.ddof = self.dim - 1
+
+        # critical chi-square value
+        self.confidence_level = confidence_level
+        self.chi2_critical_ = chi2.ppf(confidence_level, df=self.ddof)
+        self._samples = samples
+        self.alpha = 1.-confidence_level
+
+    @property
+    def samples(self):
+        return self._samples
+
+    def point_estimate(self):
         """
-        return softmax(X, axis=-1)
+        Returns the point estimate, the center of the ellipse.
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        return self.mean_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        return within_ellipse_prop(true_value, self.mean_, self.precision_matrix_, self.chi2_critical_)
+
+    def closest_point_in_region(self, p, tol=1e-6, max_iter=30):
+        return closest_point_on_ellipsoid(
+            p,
+            mean=self.mean_,
+            cov=self.cov_,
+            chi2_critical=self.chi2_critical_
+        )
+
+
+class ConfidenceEllipseTransformed(ConfidenceRegionABC):
+    """
+    Instantiates a Confidence Ellipse in a transformed space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+
+    def __init__(self, samples, transformation: CompositionalTransformation, confidence_level=0.95):
+        samples = np.asarray(samples)
+        self.transformation = transformation
+        self.confidence_level = confidence_level
+        Z = self.transformation(samples)
+        self.mean_ = np.mean(samples, axis=0)
+        # self.mean_ = self.transformation.inverse(np.mean(Z, axis=0))
+        self.conf_region_z = ConfidenceEllipseSimplex(Z, confidence_level=confidence_level)
+        self._samples = samples
+        self.alpha = 1.-confidence_level
+
+    @property
+    def samples(self):
+        return self._samples
+
+    def point_estimate(self):
+        """
+        Returns the point estimate, the center of the ellipse.
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        # The inverse of the CLR does not coincide with the true mean, because the geometric mean
+        # requires smoothing the prevalence vectors and this affects the softmax (inverse);
+        # return self.clr.inverse(self.mean_) # <- does not coincide
+        return self.mean_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        transformed_values = self.transformation(true_value)
+        return self.conf_region_z.coverage(transformed_values)
+
+    def closest_point_in_region(self, p, tol=1e-6, max_iter=30):
+        p_prime = self.transformation(p)
+        b_prime = self.conf_region_z.closest_point_in_region(p_prime, tol=tol, max_iter=max_iter)
+        b = self.transformation.inverse(b_prime)
+        return b
+
+
+class ConfidenceEllipseCLR(ConfidenceEllipseTransformed):
+    """
+    Instantiates a Confidence Ellipse in the Centered-Log Ratio (CLR) space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+    def __init__(self, samples, confidence_level=0.95):
+        super().__init__(samples, CLRtransformation(), confidence_level=confidence_level)
+
+
+class ConfidenceEllipseILR(ConfidenceEllipseTransformed):
+    """
+    Instantiates a Confidence Ellipse in the Isometric-Log Ratio (CLR) space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    """
+    def __init__(self, samples, confidence_level=0.95):
+        super().__init__(samples, ILRtransformation(), confidence_level=confidence_level)
+
+
+
+class ConfidenceIntervalsTransformed(ConfidenceRegionABC):
+    """
+    Instantiates a Confidence Interval region in a transformed space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    :param bonferroni_correction: bool (default False), if True, a Bonferroni correction
+        is applied to the significance level (`alpha`) before computing confidence intervals.
+        The correction consists of replacing `alpha` with `alpha/n_classes`. When
+        `n_classes=2` the correction is not applied because there is only one verification test
+        since the other class is constrained. This is not necessarily true for n_classes>2.
+    """
+
+    def __init__(self, samples, transformation: CompositionalTransformation, confidence_level=0.95, bonferroni_correction=False):
+        samples = np.asarray(samples)
+        self.transformation = transformation
+        self.confidence_level = confidence_level
+        Z = self.transformation(samples)
+        self.mean_ = np.mean(samples, axis=0)
+        # self.mean_ = self.transformation.inverse(np.mean(Z, axis=0))
+        self.conf_region_z = ConfidenceIntervals(Z, confidence_level=confidence_level, bonferroni_correction=bonferroni_correction)
+        self._samples = samples
+        self.alpha = 1.-confidence_level
+
+    @property
+    def samples(self):
+        return self._samples
+
+    def point_estimate(self):
+        """
+        Returns the point estimate, the center of the ellipse.
+
+        :return: np.ndarray of shape (n_classes,)
+        """
+        # The inverse of the CLR does not coincide with the true mean, because the geometric mean
+        # requires smoothing the prevalence vectors and this affects the softmax (inverse);
+        # return self.clr.inverse(self.mean_) # <- does not coincide
+        return self.mean_
+
+    def coverage(self, true_value):
+        """
+        Checks whether a value, or a sets of values, are contained in the confidence region. The method computes the
+        fraction of these that are contained in the region, if more than one value is passed. If only one value is
+        passed, then it either returns 1.0 or 0.0, for indicating the value is in the region or not, respectively.
+
+        :param true_value: a np.ndarray of shape (n_classes,) or shape (n_values, n_classes,)
+        :return: float in [0,1]
+        """
+        transformed_values = self.transformation(true_value)
+        return self.conf_region_z.coverage(transformed_values)
+
+    def coverage_soft(self, true_value):
+        transformed_values = self.transformation(true_value)
+        return self.conf_region_z.coverage_soft(transformed_values)
+
+    def winkler_scores(self, true_prev, alpha=None, add_ae=False):
+        transformed_values = self.transformation(true_prev)
+        return self.conf_region_z.winkler_scores(transformed_values, alpha=alpha, add_ae=add_ae)
+
+    def mean_winkler_score(self, true_prev, alpha=None, add_ae=False):
+        transformed_values = self.transformation(true_prev)
+        return self.conf_region_z.mean_winkler_score(transformed_values, alpha=alpha, add_ae=add_ae)
+
+
+class ConfidenceIntervalsCLR(ConfidenceIntervalsTransformed):
+    """
+    Instantiates a Confidence Intervals in the Centered-Log Ratio (CLR) space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    :param bonferroni_correction: bool (default False), if True, a Bonferroni correction
+        is applied to the significance level (`alpha`) before computing confidence intervals.
+        The correction consists of replacing `alpha` with `alpha/n_classes`. When
+        `n_classes=2` the correction is not applied because there is only one verification test
+        since the other class is constrained. This is not necessarily true for n_classes>2.
+    """
+    def __init__(self, samples, confidence_level=0.95, bonferroni_correction=False):
+        super().__init__(samples, CLRtransformation(), confidence_level=confidence_level, bonferroni_correction=bonferroni_correction)
+
+
+class ConfidenceIntervalsILR(ConfidenceIntervalsTransformed):
+    """
+    Instantiates a Confidence Intervals in the Isometric-Log Ratio (CLR) space.
+
+    :param samples: np.ndarray of shape (n_bootstrap_samples, n_classes)
+    :param confidence_level: float, the confidence level (default 0.95)
+    :param bonferroni_correction: bool (default False), if True, a Bonferroni correction
+        is applied to the significance level (`alpha`) before computing confidence intervals.
+        The correction consists of replacing `alpha` with `alpha/n_classes`. When
+        `n_classes=2` the correction is not applied because there is only one verification test
+        since the other class is constrained. This is not necessarily true for n_classes>2.
+    """
+    def __init__(self, samples, confidence_level=0.95, bonferroni_correction=False):
+        super().__init__(samples, ILRtransformation(), confidence_level=confidence_level, bonferroni_correction=bonferroni_correction)
+
 
 
 class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
@@ -399,7 +700,8 @@ class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
                  n_test_samples=500,
                  confidence_level=0.95,
                  region='intervals',
-                 random_state=None):
+                 random_state=None,
+                 verbose=False):
 
         assert isinstance(quantifier, AggregativeQuantifier), \
             f'base quantifier does not seem to be an instance of {AggregativeQuantifier.__name__}'
@@ -416,32 +718,45 @@ class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
         self.confidence_level = confidence_level
         self.region = region
         self.random_state = random_state
+        self.verbose = verbose
 
     def aggregation_fit(self, classif_predictions, labels):
-        data = LabelledCollection(classif_predictions, labels, classes=self.classes_)
+
         self.quantifiers = []
         if self.n_train_samples==1:
             self.quantifier.aggregation_fit(classif_predictions, labels)
             self.quantifiers.append(self.quantifier)
         else:
-            # model-based bootstrap (only on the aggregative part)
-            n_examples = len(data)
-            full_index = np.arange(n_examples)
-            with qp.util.temp_seed(self.random_state):
-                for i in range(self.n_train_samples):
-                    quantifier = copy.deepcopy(self.quantifier)
-                    index = resample(full_index, n_samples=n_examples)
-                    classif_predictions_i = classif_predictions.sampling_from_index(index)
-                    data_i = data.sampling_from_index(index)
-                    quantifier.aggregation_fit(classif_predictions_i, data_i)
-                    self.quantifiers.append(quantifier)
+            if classif_predictions is None or labels is None:
+                # The entire dataset was consumed for classifier training, implying there is no need for training
+                # an aggregation function. If the bootstrap method was configured to train different aggregators
+                # (i.e., self.n_train_samples>1), then an error is raise. Otherwise, the method ends.
+                if self.n_train_samples > 1:
+                    raise ValueError(
+                        f'The underlying quantifier ({self.quantifier.__class__.__name__}) has consumed, all training '
+                        f'data, meaning the aggregation function needs none, but {self.n_train_samples=} is > 1, which '
+                        f'is inconsistent.'
+                    )
+            else:
+                # model-based bootstrap (only on the aggregative part)
+                data = LabelledCollection(classif_predictions, labels, classes=self.classes_)
+                n_examples = len(data)
+                full_index = np.arange(n_examples)
+                with qp.util.temp_seed(self.random_state):
+                    for i in range(self.n_train_samples):
+                        quantifier = copy.deepcopy(self.quantifier)
+                        index = resample(full_index, n_samples=n_examples)
+                        classif_predictions_i = classif_predictions.sampling_from_index(index)
+                        data_i = data.sampling_from_index(index)
+                        quantifier.aggregation_fit(classif_predictions_i, data_i)
+                        self.quantifiers.append(quantifier)
         return self
 
     def aggregate(self, classif_predictions: np.ndarray):
         prev_mean, self.confidence = self.aggregate_conf(classif_predictions)
         return prev_mean
 
-    def aggregate_conf(self, classif_predictions: np.ndarray, confidence_level=None):
+    def aggregate_conf_sequential__(self, classif_predictions: np.ndarray, confidence_level=None):
         if confidence_level is None:
             confidence_level = self.confidence_level
 
@@ -449,7 +764,7 @@ class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
         prevs = []
         with qp.util.temp_seed(self.random_state):
             for quantifier in self.quantifiers:
-                for i in range(self.n_test_samples):
+                for i in tqdm(range(self.n_test_samples), desc='resampling', total=self.n_test_samples, disable=not self.verbose):
                     sample_i = resample(classif_predictions, n_samples=n_samples)
                     prev_i = quantifier.aggregate(sample_i)
                     prevs.append(prev_i)
@@ -459,6 +774,26 @@ class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
 
         return prev_estim, conf
 
+    def aggregate_conf(self, classif_predictions: np.ndarray, confidence_level=None):
+        confidence_level = confidence_level or self.confidence_level
+
+
+        n_samples = classif_predictions.shape[0]
+        prevs = []
+        with qp.util.temp_seed(self.random_state):
+            for quantifier in self.quantifiers:
+                results = Parallel(n_jobs=-1)(
+                    delayed(bootstrap_once)(i, classif_predictions, quantifier, n_samples)
+                    for i in range(self.n_test_samples)
+                )
+                prevs.extend(results)
+
+        prevs = np.array(prevs)
+        conf = WithConfidenceABC.construct_region(prevs, confidence_level, method=self.region)
+        prev_estim = conf.point_estimate()
+
+        return prev_estim, conf
+
     def fit(self, X, y):
         self.quantifier._check_init_parameters()
         classif_predictions, labels = self.quantifier.classifier_fit_predict(X, y)
@@ -477,6 +812,13 @@ class AggregativeBootstrap(WithConfidenceABC, AggregativeQuantifier):
         return self.quantifier._classifier_method()
 
 
+def bootstrap_once(i, classif_predictions, quantifier, n_samples):
+    idx = np.random.randint(0, len(classif_predictions), n_samples)
+    sample = classif_predictions[idx]
+    prev = quantifier.aggregate(sample)
+    return prev
+
+
 class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
     """
     `Bayesian quantification <https://arxiv.org/abs/2302.09159>`_ method (by Albert Ziegler and Paweł Czyż),
@@ -506,6 +848,8 @@ class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
     :param region: string, set to `intervals` for constructing confidence intervals (default), or to
         `ellipse` for constructing an ellipse in the probability simplex, or to `ellipse-clr` for
         constructing an ellipse in the Centered-Log Ratio (CLR) unconstrained space.
+    :param prior: an array-like with the alpha parameters of a Dirichlet prior, a scalar real value
+        to be broadcast to all classes, or the string 'uniform' for a uniform, uninformative prior (default)
     """
     def __init__(self,
                  classifier: BaseEstimator=None,
@@ -515,14 +859,21 @@ class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
                  num_samples: int = 1_000,
                  mcmc_seed: int = 0,
                  confidence_level: float = 0.95,
-                 region: str = 'intervals'):
+                 region: str = 'intervals',
+                 temperature = 1., 
+                 prior = 'uniform'):
 
         if num_warmup <= 0:
             raise ValueError(f'parameter {num_warmup=} must be a positive integer')
         if num_samples <= 0:
             raise ValueError(f'parameter {num_samples=} must be a positive integer')
+        assert ((isinstance(prior, str) and prior == 'uniform') or
+                isinstance(prior, Number) or
+                (isinstance(prior, Iterable) and all(isinstance(v, Number) for v in prior))), \
+            f'wrong type for {prior=}; expected "uniform", a real scalar, or an array-like of real values'
 
-        if _bayesian.DEPENDENCIES_INSTALLED is False:
+        bayesian = _get_bayesian_module()
+        if bayesian.DEPENDENCIES_INSTALLED is False:
             raise ImportError("Auxiliary dependencies are required. "
                               "Run `$ pip install quapy[bayes]` to install them.")
 
@@ -532,6 +883,8 @@ class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
         self.mcmc_seed = mcmc_seed
         self.confidence_level = confidence_level
         self.region = region
+        self.temperature = temperature
+        self.prior = prior
 
         # Array of shape (n_classes, n_predicted_classes,) where entry (y, c) is the number of instances
         # labeled as class y and predicted as class c.
@@ -562,11 +915,26 @@ class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
 
         n_c_unlabeled = F.counts_from_labels(classif_predictions, self.classifier.classes_).astype(float)
 
-        self._samples = _bayesian.sample_posterior(
+        n_classes = len(self.classifier.classes_)
+        if isinstance(self.prior, str) and self.prior == 'uniform':
+            alpha = np.ones(n_classes, dtype=float)
+        elif isinstance(self.prior, Number):
+            alpha = np.full(n_classes, float(self.prior), dtype=float)
+        else:
+            alpha = np.asarray(self.prior, dtype=float)
+            if alpha.ndim != 1 or len(alpha) != n_classes:
+                raise ValueError(
+                    f'wrong shape for prior; expected {n_classes} values, found shape {alpha.shape}'
+                )
+
+        bayesian = _get_bayesian_module()
+        self._samples = bayesian.sample_posterior_bayesianCC(
             n_c_unlabeled=n_c_unlabeled,
             n_y_and_c_labeled=self._n_and_c_labeled,
             num_warmup=self.num_warmup,
             num_samples=self.num_samples,
+            alpha=alpha,
+            temperature=self.temperature,
             seed=self.mcmc_seed,
         )
         return self._samples
@@ -574,15 +942,18 @@ class BayesianCC(AggregativeCrispQuantifier, WithConfidenceABC):
     def get_prevalence_samples(self):
         if self._samples is None:
             raise ValueError("sample_from_posterior must be called before get_prevalence_samples")
-        return self._samples[_bayesian.P_TEST_Y]
+        bayesian = _get_bayesian_module()
+        return self._samples[bayesian.P_TEST_Y]
 
     def get_conditional_probability_samples(self):
         if self._samples is None:
             raise ValueError("sample_from_posterior must be called before get_conditional_probability_samples")
-        return self._samples[_bayesian.P_C_COND_Y]
+        bayesian = _get_bayesian_module()
+        return self._samples[bayesian.P_C_COND_Y]
 
     def aggregate(self, classif_predictions):
-        samples = self.sample_from_posterior(classif_predictions)[_bayesian.P_TEST_Y]
+        bayesian = _get_bayesian_module()
+        samples = self.sample_from_posterior(classif_predictions)[bayesian.P_TEST_Y]
         return np.asarray(samples.mean(axis=0), dtype=float)
 
     def predict_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
@@ -637,17 +1008,19 @@ class PQ(AggregativeSoftQuantifier, BinaryAggregativeQuantifier):
         if num_samples <= 0:
             raise ValueError(f'parameter {num_samples=} must be a positive integer')
 
-        if not _bayesian.DEPENDENCIES_INSTALLED:
+        bayesian = _get_bayesian_module()
+        if not bayesian.DEPENDENCIES_INSTALLED:
             raise ImportError("Auxiliary dependencies are required. "
                               "Run `$ pip install quapy[bayes]` to install them.")
 
         super().__init__(classifier, fit_classifier, val_split)
+        
         self.nbins = nbins
         self.fixed_bins = fixed_bins
         self.num_warmup = num_warmup
         self.num_samples = num_samples
         self.stan_seed = stan_seed
-        self.stan_code = _bayesian.load_stan_file()
+        self.stan_code = bayesian.load_stan_file()
         self.confidence_level = confidence_level
         self.region = region
 
@@ -676,7 +1049,8 @@ class PQ(AggregativeSoftQuantifier, BinaryAggregativeQuantifier):
     def aggregate(self, classif_predictions):
         Px_test = classif_predictions[:, self.pos_label]
         test_hist, _ = np.histogram(Px_test, bins=self.bin_limits)
-        prevs = _bayesian.pq_stan(
+        bayesian = _get_bayesian_module()
+        prevs = bayesian.pq_stan(
             self.stan_code, self.nbins, self.pos_hist, self.neg_hist, test_hist,
             self.num_samples, self.num_warmup, self.stan_seed
         ).flatten()
@@ -694,5 +1068,3 @@ class PQ(AggregativeSoftQuantifier, BinaryAggregativeQuantifier):
     def predict_conf(self, instances, confidence_level=None) -> (np.ndarray, ConfidenceRegionABC):
         predictions = self.classify(instances)
         return self.aggregate_conf(predictions, confidence_level=confidence_level)
-
-

quapy/protocol.py

@@ -10,6 +10,8 @@ from quapy.data import LabelledCollection
 import quapy.functional as F
 from os.path import exists
 from glob import glob
+from collections.abc import Iterable
+from numbers import Number
 
 
 class AbstractProtocol(metaclass=ABCMeta):
@@ -171,7 +173,7 @@ class AbstractStochasticSeededProtocol(AbstractProtocol):
         return sample
 
 
-class OnLabelledCollectionProtocol:
+class OnLabelledCollectionProtocol(AbstractStochasticSeededProtocol):
     """
     Protocols that generate samples from a :class:`qp.data.LabelledCollection` object.
     """
@@ -229,8 +231,17 @@ class OnLabelledCollectionProtocol:
         elif return_type=='index':
             return lambda lc,params:params
 
+    def sample(self, index):
+        """
+        Realizes the sample given the index of the instances.
 
-class APP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
+        :param index: indexes of the instances to select
+        :return: an instance of :class:`qp.data.LabelledCollection`
+        """
+        return self.data.sampling_from_index(index)
+
+
+class APP(OnLabelledCollectionProtocol):
     """
     Implementation of the artificial prevalence protocol (APP).
     The APP consists of exploring a grid of prevalence values containing `n_prevalences` points (e.g.,
@@ -311,15 +322,6 @@ class APP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
             indexes.append(index)
         return indexes
 
-    def sample(self, index):
-        """
-        Realizes the sample given the index of the instances.
-
-        :param index: indexes of the instances to select
-        :return: an instance of :class:`qp.data.LabelledCollection`
-        """
-        return self.data.sampling_from_index(index)
-
     def total(self):
         """
         Returns the number of samples that will be generated
@@ -329,7 +331,7 @@ class APP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
         return F.num_prevalence_combinations(self.n_prevalences, self.data.n_classes, self.repeats)
 
 
-class NPP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
+class NPP(OnLabelledCollectionProtocol):
     """
     A generator of samples that implements the natural prevalence protocol (NPP). The NPP consists of drawing
     samples uniformly at random, therefore approximately preserving the natural prevalence of the collection.
@@ -365,15 +367,6 @@ class NPP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
             indexes.append(index)
         return indexes
 
-    def sample(self, index):
-        """
-        Realizes the sample given the index of the instances.
-
-        :param index: indexes of the instances to select
-        :return: an instance of :class:`qp.data.LabelledCollection`
-        """
-        return self.data.sampling_from_index(index)
-
     def total(self):
         """
         Returns the number of samples that will be generated (equals to "repeats")
@@ -383,7 +376,7 @@ class NPP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
         return self.repeats
 
 
-class UPP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
+class UPP(OnLabelledCollectionProtocol):
     """
     A variant of :class:`APP` that, instead of using a grid of equidistant prevalence values,
     relies on the Kraemer algorithm for sampling unit (k-1)-simplex uniformly at random, with
@@ -423,14 +416,69 @@ class UPP(AbstractStochasticSeededProtocol, OnLabelledCollectionProtocol):
             indexes.append(index)
         return indexes
 
-    def sample(self, index):
+    def total(self):
         """
-        Realizes the sample given the index of the instances.
+        Returns the number of samples that will be generated (equals to "repeats")
 
-        :param index: indexes of the instances to select
-        :return: an instance of :class:`qp.data.LabelledCollection`
+        :return: int
         """
-        return self.data.sampling_from_index(index)
+        return self.repeats
+
+
+class DirichletProtocol(OnLabelledCollectionProtocol):
+    """
+    A protocol that establishes a prior Dirichlet distribution for the prevalence of the samples.
+    Note that providing an all-ones vector of Dirichlet parameters is equivalent to invoking the
+    APP protocol (although each protocol will generate a different series of samples given a
+    fixed seed, since the implementation is different).
+
+    :param data: a `LabelledCollection` from which the samples will be drawn
+    :param alpha: an array-like of shape (n_classes,) with the parameters of the Dirichlet distribution
+    :param sample_size: integer, the number of instances in each sample; if None (default) then it is taken from
+        qp.environ["SAMPLE_SIZE"]. If this is not set, a ValueError exception is raised.
+    :param repeats: the number of samples to generate. Default is 100.
+    :param random_state: allows replicating samples across runs (default 0, meaning that the sequence of samples
+        will be the same every time the protocol is called)
+    :param return_type: set to "sample_prev" (default) to get the pairs of (sample, prevalence) at each iteration, or
+        to "labelled_collection" to get instead instances of LabelledCollection
+    """
+
+    def __init__(self, data: LabelledCollection, alpha, sample_size=None, repeats=100, random_state=0,
+                 return_type='sample_prev'):
+        #assert ((isinstance(alpha, str) and alpha == 'uniform') or
+        #        isinstance(alpha, Number) or
+        #        (isinstance(alpha, Iterable) and all(isinstance(v, Number) for v in alpha))), \
+        #    f'wrong type for {alpha=}; expected "uniform", a real scalar, or an array-like of real values'
+        
+        n_classes = data.n_classes
+        if isinstance(alpha, str) and alpha == 'uniform':
+            self.alpha = np.ones(n_classes, dtype=float)
+        elif isinstance(alpha, Number):
+            self.alpha = np.full(n_classes, float(alpha), dtype=float)
+        else:
+            self.alpha = np.asarray(alpha, dtype=float)
+            if self.alpha.ndim != 1 or len(self.alpha) != n_classes:
+                raise ValueError(
+                    f'wrong shape for alpha; expected {n_classes} values, found shape {self.alpha.shape}'
+                )
+
+        super(DirichletProtocol, self).__init__(random_state)
+        self.data = data
+        #self.alpha = alpha
+        self.sample_size = qp._get_sample_size(sample_size)
+        self.repeats = repeats
+        self.random_state = random_state
+        self.collator = OnLabelledCollectionProtocol.get_collator(return_type)
+
+    def samples_parameters(self):
+        """
+        Return all the necessary parameters to replicate the samples.
+
+        :return: a list of indexes that realize the sampling
+        """
+        prevs = np.random.dirichlet(self.alpha, size=self.repeats)
+        indexes = [self.data.sampling_index(self.sample_size, *prevs_i) for prevs_i in prevs]
+        return indexes
 
     def total(self):
         """
@@ -450,7 +498,7 @@ class DomainMixer(AbstractStochasticSeededProtocol):
     :param sample_size: integer, the number of instances in each sample; if None (default) then it is taken from
         qp.environ["SAMPLE_SIZE"]. If this is not set, a ValueError exception is raised.
     :param repeats: int, number of samples to draw for every mixture rate
-    :param prevalence: the prevalence to preserv along the mixtures. If specified, should be an array containing
+    :param prevalence: the prevalence to preserve along the mixtures. If specified, should be an array containing
         one prevalence value (positive float) for each class and summing up to one. If not specified, the prevalence
         will be taken from the domain A (default).
     :param mixture_points: an integer indicating the number of points to take from a linear scale (e.g., 21 will

quapy/tests/__init__.py

@@ -0,0 +1,10 @@
+"""
+Fast unit tests live in ``test_*.py`` and are intended to run without network
+access or large external resources.
+
+Slow integration tests that download datasets or depend on optional stacks live
+in ``integration_*.py`` and are meant to be run explicitly, e.g.:
+
+    python -m unittest quapy.tests.integration_datasets
+    python -m unittest quapy.tests.integration_methods
+"""

quapy/tests/_synthetic.py

@@ -0,0 +1,48 @@
+import numpy as np
+from sklearn.datasets import make_classification
+
+from quapy.data import LabelledCollection
+from quapy.data.base import Dataset
+
+
+def make_labelled_collection(
+    n_samples=200,
+    n_features=12,
+    n_classes=2,
+    class_sep=1.5,
+    random_state=0,
+):
+    n_informative = min(n_features, max(4, n_classes * 2))
+    X, y = make_classification(
+        n_samples=n_samples,
+        n_features=n_features,
+        n_informative=n_informative,
+        n_redundant=0,
+        n_repeated=0,
+        n_classes=n_classes,
+        n_clusters_per_class=1,
+        class_sep=class_sep,
+        random_state=random_state,
+    )
+    classes = np.arange(n_classes)
+    return LabelledCollection(X, y, classes=classes)
+
+
+def make_dataset(
+    n_train=150,
+    n_test=80,
+    n_features=12,
+    n_classes=2,
+    class_sep=1.5,
+    random_state=0,
+    name='synthetic',
+):
+    data = make_labelled_collection(
+        n_samples=n_train + n_test,
+        n_features=n_features,
+        n_classes=n_classes,
+        class_sep=class_sep,
+        random_state=random_state,
+    )
+    training, test = data.split_stratified(train_prop=n_train / (n_train + n_test), random_state=random_state)
+    return Dataset(training, test, name=name)

quapy/tests/integration_datasets.py

@@ -1,3 +1,10 @@
+"""
+Integration tests for dataset fetchers and large external resources.
+
+This module is intentionally excluded from default ``unittest`` discovery by
+using an ``integration_*.py`` filename.
+"""
+
 import os
 import unittest
 
@@ -5,11 +12,11 @@ from sklearn.feature_extraction.text import TfidfVectorizer
 from sklearn.linear_model import LogisticRegression
 
 import quapy.functional as F
-from quapy.method.aggregative import PCC
 from quapy.data.datasets import *
+from quapy.method.aggregative import PCC
 
 
-class TestDatasets(unittest.TestCase):
+class IntegrationDatasetsTest(unittest.TestCase):
 
     def new_quantifier(self):
         return PCC(LogisticRegression(C=0.001, max_iter=100))
@@ -17,13 +24,11 @@ class TestDatasets(unittest.TestCase):
     def _check_dataset(self, dataset):
         train, test = dataset.reduce().train_test
         q = self.new_quantifier()
-        print(f'testing method {q} in {dataset.name}...', end='')
-        if len(train)>500:
+        if len(train) > 500:
             train = train.sampling(500)
         q.fit(*dataset.training.Xy)
         estim_prevalences = q.predict(dataset.test.instances)
         self.assertTrue(F.check_prevalence_vector(estim_prevalences))
-        print(f'[done]')
 
     def _check_samples(self, gen, q, max_samples_test=5, vectorizer=None):
         for X, p in gen():
@@ -37,54 +42,37 @@ class TestDatasets(unittest.TestCase):
 
     def test_reviews(self):
         for dataset_name in REVIEWS_SENTIMENT_DATASETS:
-            print(f'loading dataset {dataset_name}...', end='')
             dataset = fetch_reviews(dataset_name, tfidf=True, min_df=10)
-            dataset.stats()
             dataset.reduce()
-            print(f'[done]')
             self._check_dataset(dataset)
 
     def test_twitter(self):
-        # all the datasets are contained in the same resource; if the first one
-        # works, there is no need to test for the rest
         for dataset_name in TWITTER_SENTIMENT_DATASETS_TEST[:1]:
-            print(f'loading dataset {dataset_name}...', end='')
             dataset = fetch_twitter(dataset_name, min_df=10)
-            dataset.stats()
             dataset.reduce()
-            print(f'[done]')
             self._check_dataset(dataset)
 
     def test_UCIBinaryDataset(self):
         for dataset_name in UCI_BINARY_DATASETS:
-            print(f'loading dataset {dataset_name}...', end='')
             dataset = fetch_UCIBinaryDataset(dataset_name)
-            dataset.stats()
             dataset.reduce()
-            print(f'[done]')
             self._check_dataset(dataset)
 
     def test_UCIMultiDataset(self):
         for dataset_name in UCI_MULTICLASS_DATASETS:
-            print(f'loading dataset {dataset_name}...', end='')
             dataset = fetch_UCIMulticlassDataset(dataset_name)
-            dataset.stats()
             n_classes = dataset.n_classes
             uniform_prev = F.uniform_prevalence(n_classes)
             dataset.training = dataset.training.sampling(100, *uniform_prev)
             dataset.test = dataset.test.sampling(100, *uniform_prev)
-            print(f'[done]')
             self._check_dataset(dataset)
 
     def test_lequa2022(self):
         if os.environ.get('QUAPY_TESTS_OMIT_LARGE_DATASETS'):
-            print("omitting test_lequa2022 because QUAPY_TESTS_OMIT_LARGE_DATASETS is set")
             return
 
         for dataset_name in LEQUA2022_VECTOR_TASKS:
-            print(f'LeQu2022: loading dataset {dataset_name}...', end='')
             train, gen_val, gen_test = fetch_lequa2022(dataset_name)
-            train.stats()
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             q = self.new_quantifier()
@@ -93,9 +81,7 @@ class TestDatasets(unittest.TestCase):
             self._check_samples(gen_test, q, max_samples_test=5)
 
         for dataset_name in LEQUA2022_TEXT_TASKS:
-            print(f'LeQu2022: loading dataset {dataset_name}...', end='')
             train, gen_val, gen_test = fetch_lequa2022(dataset_name)
-            train.stats()
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             tfidf = TfidfVectorizer()
@@ -107,13 +93,10 @@ class TestDatasets(unittest.TestCase):
 
     def test_lequa2024(self):
         if os.environ.get('QUAPY_TESTS_OMIT_LARGE_DATASETS'):
-            print("omitting test_lequa2024 because QUAPY_TESTS_OMIT_LARGE_DATASETS is set")
             return
 
         for task in LEQUA2024_TASKS:
-            print(f'LeQu2024: loading task {task}...', end='')
             train, gen_val, gen_test = fetch_lequa2024(task, merge_T3=True)
-            train.stats()
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             q = self.new_quantifier()
@@ -121,16 +104,12 @@ class TestDatasets(unittest.TestCase):
             self._check_samples(gen_val, q, max_samples_test=5)
             self._check_samples(gen_test, q, max_samples_test=5)
 
-
     def test_IFCB(self):
         if os.environ.get('QUAPY_TESTS_OMIT_LARGE_DATASETS'):
-            print("omitting test_IFCB because QUAPY_TESTS_OMIT_LARGE_DATASETS is set")
             return
 
-        print(f'loading dataset IFCB.')
         for mod_sel in [False, True]:
             train, gen = fetch_IFCB(single_sample_train=True, for_model_selection=mod_sel)
-            train.stats()
             n_classes = train.n_classes
             train = train.sampling(100, *F.uniform_prevalence(n_classes))
             q = self.new_quantifier()

quapy/tests/integration_methods.py

@@ -0,0 +1,42 @@
+"""
+Integration tests for optional or resource-heavy method end-to-end checks.
+
+This module is intentionally excluded from default ``unittest`` discovery by
+using an ``integration_*.py`` filename.
+"""
+
+import unittest
+
+import quapy as qp
+from quapy.functional import check_prevalence_vector
+
+
+class IntegrationMethodsTest(unittest.TestCase):
+
+    def test_quanet(self):
+        try:
+            import quapy.classification.neural
+        except ModuleNotFoundError:
+            print('the torch package is not installed; skipping integration test for QuaNet')
+            return
+
+        qp.environ['SAMPLE_SIZE'] = 10
+
+        dataset = qp.datasets.fetch_reviews('kindle', pickle=True).reduce()
+        qp.data.preprocessing.index(dataset, min_df=5, inplace=True)
+
+        from quapy.classification.neural import CNNnet
+        from quapy.classification.neural import NeuralClassifierTrainer
+        from quapy.method.meta import QuaNet
+
+        cnn = CNNnet(dataset.vocabulary_size, dataset.n_classes)
+        learner = NeuralClassifierTrainer(cnn, device='cpu')
+        model = QuaNet(learner, device='cpu', n_epochs=2, tr_iter_per_poch=10, va_iter_per_poch=10, patience=2)
+
+        model.fit(*dataset.training.Xy)
+        estim_prevalences = model.predict(dataset.test.instances)
+        self.assertTrue(check_prevalence_vector(estim_prevalences))
+
+
+if __name__ == '__main__':
+    unittest.main()

quapy/tests/test_evaluation.py

@@ -1,43 +1,41 @@
+import inspect
 import unittest
-
-import numpy as np
-
-import quapy as qp
-from sklearn.linear_model import LogisticRegression
 from time import time
 
+import numpy as np
+from sklearn.linear_model import LogisticRegression
+
+import quapy as qp
 from quapy.error import QUANTIFICATION_ERROR_SINGLE_NAMES
 from quapy.method.aggregative import EMQ, PCC
 from quapy.method.base import BaseQuantifier
+from quapy.tests._synthetic import make_dataset
 
 
 class EvalTestCase(unittest.TestCase):
 
+    @classmethod
+    def setUpClass(cls):
+        cls.data = make_dataset(n_train=140, n_test=90, n_classes=2, random_state=7, name='eval')
+
     def test_eval_speedup(self):
-        """
-        Checks whether the speed-up heuristics used by qp.evaluation work, i.e., actually save time
-        """
-
-        data = qp.datasets.fetch_reviews('hp', tfidf=True, min_df=10, pickle=True)
-        train, test = data.training, data.test
-
-        protocol = qp.protocol.APP(test, sample_size=1000, n_prevalences=11, repeats=1, random_state=1)
+        train, test = self.data.training, self.data.test
+        protocol = qp.protocol.APP(test, sample_size=30, n_prevalences=5, repeats=1, random_state=1)
 
         class SlowLR(LogisticRegression):
             def predict_proba(self, X):
-                import time
-                time.sleep(1)
+                import time as _time
+                _time.sleep(0.05)
                 return super().predict_proba(X)
 
-        emq = EMQ(SlowLR()).fit(*train.Xy)
+        emq = EMQ(SlowLR(max_iter=1000)).fit(*train.Xy)
 
         tinit = time()
-        score = qp.evaluation.evaluate(emq, protocol, error_metric='mae', verbose=True, aggr_speedup='force')
-        tend_optim = time()-tinit
-        print(f'evaluation (with optimization) took {tend_optim}s [MAE={score:.4f}]')
+        score = qp.evaluation.evaluate(emq, protocol, error_metric='mae', aggr_speedup='force')
+        tend_optim = time() - tinit
+        self.assertTrue(isinstance(score, float))
 
         class NonAggregativeEMQ(BaseQuantifier):
-
             def __init__(self, cls):
                 self.emq = EMQ(cls)
 
@@ -48,31 +46,32 @@ class EvalTestCase(unittest.TestCase):
                 self.emq.fit(X, y)
                 return self
 
-        emq = NonAggregativeEMQ(SlowLR()).fit(*train.Xy)
+        emq = NonAggregativeEMQ(SlowLR(max_iter=1000)).fit(*train.Xy)
 
         tinit = time()
-        score = qp.evaluation.evaluate(emq, protocol, error_metric='mae', verbose=True)
+        score = qp.evaluation.evaluate(emq, protocol, error_metric='mae')
         tend_no_optim = time() - tinit
-        print(f'evaluation (w/o optimization) took {tend_no_optim}s [MAE={score:.4f}]')
-
-        self.assertEqual(tend_no_optim>(tend_optim/2), True)
+        self.assertTrue(isinstance(score, float))
+        self.assertGreater(tend_no_optim, tend_optim)
 
     def test_evaluation_output(self):
-        """
-        Checks the evaluation functions return correct types for different error_metrics
-        """
+        train, test = self.data.training, self.data.test
+        qp.environ['SAMPLE_SIZE'] = 30
+        protocol = qp.protocol.APP(test, sample_size=30, n_prevalences=5, repeats=1, random_state=0)
+        q = PCC(LogisticRegression(max_iter=1000)).fit(*train.Xy)
 
-        data = qp.datasets.fetch_reviews('hp', tfidf=True, min_df=10, pickle=True).reduce(n_train=100, n_test=100)
-        train, test = data.training, data.test
+        def supports_evaluation(err):
+            required = [
+                p for p in inspect.signature(err).parameters.values()
+                if p.kind in (p.POSITIONAL_ONLY, p.POSITIONAL_OR_KEYWORD) and p.default is inspect._empty
+            ]
+            return len(required) <= 2
 
-        qp.environ['SAMPLE_SIZE']=100
-
-        protocol = qp.protocol.APP(test, random_state=0)
-
-        q = PCC(LogisticRegression()).fit(*train.Xy)
-
-        single_errors = list(QUANTIFICATION_ERROR_SINGLE_NAMES)
-        averaged_errors = ['m'+e for e in single_errors]
+        single_errors = [
+            e for e in QUANTIFICATION_ERROR_SINGLE_NAMES
+            if supports_evaluation(qp.error.from_name(e))
+        ]
+        averaged_errors = ['m' + e for e in single_errors]
         single_errors = single_errors + [qp.error.from_name(e) for e in single_errors]
         averaged_errors = averaged_errors + [qp.error.from_name(e) for e in averaged_errors]
         for error_metric, averaged_error_metric in zip(single_errors, averaged_errors):
@@ -81,7 +80,6 @@ class EvalTestCase(unittest.TestCase):
 
             scores = qp.evaluation.evaluate(q, protocol, error_metric=error_metric)
             self.assertTrue(isinstance(scores, np.ndarray))
-
             self.assertEqual(scores.mean(), score)
 
 

quapy/tests/test_methods.py

@@ -1,125 +1,122 @@
 import itertools
+import inspect
 import unittest
 
 from sklearn.linear_model import LogisticRegression
 
-import quapy as qp
+from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_METHODS
 from quapy.method.aggregative import ACC
 from quapy.method.meta import Ensemble
-from quapy.method import AGGREGATIVE_METHODS, BINARY_METHODS, NON_AGGREGATIVE_METHODS
 from quapy.functional import check_prevalence_vector
+from quapy.tests._synthetic import make_dataset
+import quapy as qp
 
-# a random selection of composed methods to test the qunfold integration
-from quapy.method.composable import check_compatible_qunfold_version
+OPTIONAL_AGGREGATIVE_METHODS = {
+    'BayesianCC',
+    'BayesianKDEy',
+    'BayesianMAPLS',
+    'PQ',
+}
 
-from quapy.method.composable import (
-    ComposableQuantifier,
-    LeastSquaresLoss,
-    HellingerSurrogateLoss,
-    ClassRepresentation,
-    HistogramRepresentation,
-    CVClassifier
-)
-
-COMPOSABLE_METHODS = [
-    ComposableQuantifier( # ACC
-        LeastSquaresLoss(),
-        ClassRepresentation(CVClassifier(LogisticRegression()))
-    ),
-    ComposableQuantifier( # HDy
-        HellingerSurrogateLoss(),
-        HistogramRepresentation(
-            3, # 3 bins per class
-            preprocessor = ClassRepresentation(CVClassifier(LogisticRegression()))
-        )
-    ),
-]
 
 class TestMethods(unittest.TestCase):
 
-    tiny_dataset_multiclass = qp.datasets.fetch_UCIMulticlassDataset('academic-success').reduce(n_test=10)
-    tiny_dataset_binary = qp.datasets.fetch_UCIBinaryDataset('ionosphere').reduce(n_test=10)
+    tiny_dataset_multiclass = make_dataset(
+        n_train=140, n_test=40, n_classes=3, n_features=12, random_state=11, name='synthetic-multiclass'
+    )
+    tiny_dataset_binary = make_dataset(
+        n_train=140, n_test=40, n_classes=2, n_features=12, random_state=13, name='synthetic-binary'
+    )
     datasets = [tiny_dataset_binary, tiny_dataset_multiclass]
 
     def test_aggregative(self):
         for dataset in TestMethods.datasets:
-            learner = LogisticRegression()
+            learner = LogisticRegression(max_iter=2000)
             learner.fit(*dataset.training.Xy)
 
             for model in AGGREGATIVE_METHODS:
+                if model.__name__ in OPTIONAL_AGGREGATIVE_METHODS:
+                    continue
                 if not dataset.binary and model in BINARY_METHODS:
-                    print(f'skipping the test of binary model {model.__name__} on multiclass dataset {dataset.name}')
                     continue
 
-                q = model(learner, fit_classifier=False)
-                print('testing', q)
+                kwargs = {'fit_classifier': False}
+                if 'val_split' in inspect.signature(model.__init__).parameters:
+                    kwargs['val_split'] = None
+                q = model(learner, **kwargs)
                 q.fit(*dataset.training.Xy)
                 estim_prevalences = q.predict(dataset.test.X)
                 self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_non_aggregative(self):
         for dataset in TestMethods.datasets:
-
             for model in NON_AGGREGATIVE_METHODS:
                 if not dataset.binary and model in BINARY_METHODS:
-                    print(f'skipping the test of binary model {model.__name__} on multiclass dataset {dataset.name}')
                     continue
 
                 q = model()
-                print(f'testing {q} on dataset {dataset.name}')
                 q.fit(*dataset.training.Xy)
                 estim_prevalences = q.predict(dataset.test.X)
                 self.assertTrue(check_prevalence_vector(estim_prevalences))
 
     def test_ensembles(self):
-        qp.environ['SAMPLE_SIZE'] = 10
+        qp.environ['SAMPLE_SIZE'] = 20
 
-        base_quantifier = ACC(LogisticRegression())
+        def policy_supported(policy):
+            if policy in {'ave', 'ptr', 'ds'}:
+                return True
+            err = qp.error.from_name(policy)
+            required = [
+                p for p in inspect.signature(err).parameters.values()
+                if p.kind in (p.POSITIONAL_ONLY, p.POSITIONAL_OR_KEYWORD) and p.default is inspect._empty
+            ]
+            return len(required) <= 2
+
+        base_quantifier = ACC(LogisticRegression(max_iter=2000))
         for dataset, policy in itertools.product(TestMethods.datasets, Ensemble.VALID_POLICIES):
+            if not policy_supported(policy):
+                continue
             if not dataset.binary and policy == 'ds':
-                print(f'skipping the test of binary policy ds on non-binary dataset {dataset}')
                 continue
 
-            print(f'testing {base_quantifier} on dataset {dataset.name} with {policy=}')
-            ensemble = Ensemble(quantifier=base_quantifier, size=3, policy=policy, n_jobs=-1)
+            ensemble = Ensemble(quantifier=base_quantifier, size=3, policy=policy, n_jobs=1)
             ensemble.fit(*dataset.training.Xy)
             estim_prevalences = ensemble.predict(dataset.test.instances)
             self.assertTrue(check_prevalence_vector(estim_prevalences))
 
-    def test_quanet(self):
+    def test_composable(self):
         try:
-            import quapy.classification.neural
-        except ModuleNotFoundError:
-            print('the torch package is not installed; skipping unit test for QuaNet')
+            from quapy.method.composable import check_compatible_qunfold_version
+            from quapy.method.composable import (
+                ComposableQuantifier,
+                LeastSquaresLoss,
+                HellingerSurrogateLoss,
+                ClassRepresentation,
+                HistogramRepresentation,
+                CVClassifier,
+            )
+        except ImportError:
             return
 
-        qp.environ['SAMPLE_SIZE'] = 10
+        composable_methods = [
+            ComposableQuantifier(
+                LeastSquaresLoss(),
+                ClassRepresentation(CVClassifier(LogisticRegression()))
+            ),
+            ComposableQuantifier(
+                HellingerSurrogateLoss(),
+                HistogramRepresentation(
+                    3,
+                    preprocessor=ClassRepresentation(CVClassifier(LogisticRegression()))
+                )
+            ),
+        ]
 
-        # load the kindle dataset as text, and convert words to numerical indexes
-        dataset = qp.datasets.fetch_reviews('kindle', pickle=True).reduce()
-        qp.data.preprocessing.index(dataset, min_df=5, inplace=True)
-
-        from quapy.classification.neural import CNNnet
-        cnn = CNNnet(dataset.vocabulary_size, dataset.n_classes)
-
-        from quapy.classification.neural import NeuralClassifierTrainer
-        learner = NeuralClassifierTrainer(cnn, device='cpu')
-
-        from quapy.method.meta import QuaNet
-        model = QuaNet(learner, device='cpu', n_epochs=2, tr_iter_per_poch=10, va_iter_per_poch=10, patience=2)
-
-        model.fit(*dataset.training.Xy)
-        estim_prevalences = model.predict(dataset.test.instances)
-        self.assertTrue(check_prevalence_vector(estim_prevalences))
-
-    def test_composable(self):
         if check_compatible_qunfold_version():
             for dataset in TestMethods.datasets:
-                for q in COMPOSABLE_METHODS:
-                    print('testing', q)
+                for q in composable_methods:
                     q.fit(*dataset.training.Xy)
                     estim_prevalences = q.predict(dataset.test.X)
-                    print(estim_prevalences)
                     self.assertTrue(check_prevalence_vector(estim_prevalences))
         else:
             from quapy.method.composable import __old_version_message

quapy/tests/test_modsel.py

@@ -1,104 +1,87 @@
+import time
 import unittest
 
 import numpy as np
 from sklearn.linear_model import LogisticRegression
 
-import quapy as qp
 from quapy.method.aggregative import PACC
 from quapy.model_selection import GridSearchQ
 from quapy.protocol import APP
-import time
+from quapy.tests._synthetic import make_dataset
 
 
 class ModselTestCase(unittest.TestCase):
 
+    @classmethod
+    def setUpClass(cls):
+        data = make_dataset(
+            n_train=220,
+            n_test=120,
+            n_classes=2,
+            n_features=16,
+            class_sep=1.8,
+            random_state=1,
+            name='modsel',
+        )
+        cls.training, cls.validation = data.training.split_stratified(0.7, random_state=1)
+
     def test_modsel(self):
         """
-        Checks whether a model selection exploration takes a good hyperparameter
+        Checks whether a model selection exploration picks the better hyperparameter.
         """
-
         q = PACC(LogisticRegression(random_state=1, max_iter=5000))
-
-        data = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=10).reduce(random_state=1)
-        training, validation = data.training.split_stratified(0.7, random_state=1)
-
-        param_grid = {'classifier__C': [0.000001, 10.]}
-        app = APP(validation, sample_size=100, random_state=1)
+        param_grid = {'classifier__C': [0.000001, 10.0]}
+        app = APP(self.validation, sample_size=30, n_prevalences=5, repeats=1, random_state=1)
         q = GridSearchQ(
-            q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, verbose=True, n_jobs=-1
-        ).fit(*training.Xy)
-        print('best params', q.best_params_)
-        print('best score', q.best_score_)
+            q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, verbose=False, n_jobs=-1
+        ).fit(*self.training.Xy)
 
         self.assertEqual(q.best_params_['classifier__C'], 10.0)
         self.assertEqual(q.best_model().get_params()['classifier__C'], 10.0)
 
     def test_modsel_parallel(self):
         """
-        Checks whether a parallelized model selection actually is faster than a sequential exploration but
-        obtains the same optimal parameters
+        Checks whether sequential and parallel model selection agree on the best parameters.
         """
-
         q = PACC(LogisticRegression(random_state=1, max_iter=3000))
-
-        data = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=50)
-        training, validation = data.training.split_stratified(0.7, random_state=1)
-
-        param_grid = {'classifier__C': np.logspace(-3,3,7), 'classifier__class_weight': ['balanced', None]}
-        app = APP(validation, sample_size=100, random_state=1)
+        param_grid = {'classifier__C': np.logspace(-3, 3, 7), 'classifier__class_weight': ['balanced', None]}
+        app = APP(self.validation, sample_size=30, n_prevalences=5, repeats=1, random_state=1)
 
         def do_gridsearch(n_jobs):
-            print('starting model selection in sequential exploration')
             t_init = time.time()
             modsel = GridSearchQ(
-                q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, n_jobs=n_jobs, verbose=True
-            ).fit(*training.Xy)
-            t_end = time.time()-t_init
-            best_c = modsel.best_params_['classifier__C']
-            print(f'[done] took {t_end:.2f}s best C = {best_c}')
-            return t_end, best_c
+                q, param_grid, protocol=app, error='mae', refit=False, timeout=-1, n_jobs=n_jobs, verbose=False
+            ).fit(*self.training.Xy)
+            t_end = time.time() - t_init
+            return t_end, modsel.best_params_
 
-        tend_seq, best_c_seq = do_gridsearch(n_jobs=1)
-        tend_par, best_c_par = do_gridsearch(n_jobs=-1)
-
-        print(tend_seq, best_c_seq)
-        print(tend_par, best_c_par)
-
-        self.assertEqual(best_c_seq, best_c_par)
-        self.assertLess(tend_par, tend_seq)
+        _, best_seq = do_gridsearch(n_jobs=1)
+        _, best_par = do_gridsearch(n_jobs=-1)
 
+        self.assertEqual(best_seq, best_par)
 
     def test_modsel_timeout(self):
 
         class SlowLR(LogisticRegression):
             def fit(self, X, y, sample_weight=None):
-                import time
-                time.sleep(10)
-                super(SlowLR, self).fit(X, y, sample_weight)
+                time.sleep(2)
+                return super().fit(X, y, sample_weight)
 
-        q = PACC(SlowLR())
+        q = PACC(SlowLR(max_iter=1000))
+        param_grid = {'classifier__C': np.logspace(-1, 1, 3)}
+        app = APP(self.validation, sample_size=30, n_prevalences=5, repeats=1, random_state=1)
 
-        data = qp.datasets.fetch_reviews('imdb', tfidf=True, min_df=10).reduce(random_state=1)
-        training, validation = data.training.split_stratified(0.7, random_state=1)
-
-        param_grid = {'classifier__C': np.logspace(-1,1,3)}
-        app = APP(validation, sample_size=100, random_state=1)
-
-        print('Expecting TimeoutError to be raised')
         modsel = GridSearchQ(
-            q, param_grid, protocol=app, timeout=3, n_jobs=-1, verbose=True, raise_errors=True
+            q, param_grid, protocol=app, timeout=1, n_jobs=-1, verbose=False, raise_errors=True
         )
         with self.assertRaises(TimeoutError):
-            modsel.fit(*training.Xy)
+            modsel.fit(*self.training.Xy)
 
-        print('Expecting ValueError to be raised')
         modsel = GridSearchQ(
-            q, param_grid, protocol=app, timeout=3, n_jobs=-1, verbose=True, raise_errors=False
+            q, param_grid, protocol=app, timeout=1, n_jobs=-1, verbose=False, raise_errors=False
         )
         with self.assertRaises(ValueError):
-            # this exception is not raised because of the timeout, but because no combination of hyperparams
-            # succedded (in this case, a ValueError is raised, regardless of "raise_errors"
-            modsel.fit(*training.Xy)
+            modsel.fit(*self.training.Xy)
 
 
 if __name__ == '__main__':

quapy/tests/test_protocols.py

@@ -3,7 +3,7 @@ import numpy as np
 
 import quapy.functional
 from quapy.data import LabelledCollection
-from quapy.protocol import APP, NPP, UPP, DomainMixer, AbstractStochasticSeededProtocol
+from quapy.protocol import APP, NPP, UPP, DomainMixer, AbstractStochasticSeededProtocol, DirichletProtocol
 
 
 def mock_labelled_collection(prefix=''):
@@ -138,6 +138,31 @@ class TestProtocols(unittest.TestCase):
 
         self.assertNotEqual(samples1, samples2)
 
+    def test_dirichlet_replicate(self):
+        data = mock_labelled_collection()
+        p = DirichletProtocol(data, alpha=[1, 2, 3, 4], sample_size=5, repeats=10, random_state=42)
+
+        samples1 = samples_to_str(p)
+        samples2 = samples_to_str(p)
+
+        self.assertEqual(samples1, samples2)
+
+        p = DirichletProtocol(data, alpha=[1, 2, 3, 4], sample_size=5, repeats=10, random_state=0)
+
+        samples1 = samples_to_str(p)
+        samples2 = samples_to_str(p)
+
+        self.assertEqual(samples1, samples2)
+
+    def test_dirichlet_not_replicate(self):
+        data = mock_labelled_collection()
+        p = DirichletProtocol(data, alpha=[1, 2, 3, 4], sample_size=5, repeats=10, random_state=None)
+
+        samples1 = samples_to_str(p)
+        samples2 = samples_to_str(p)
+
+        self.assertNotEqual(samples1, samples2)
+
     def test_covariate_shift_replicate(self):
         dataA = mock_labelled_collection('domA')
         dataB = mock_labelled_collection('domB')

quapy/tests/test_replicability.py

@@ -1,19 +1,29 @@
 import unittest
+
+import numpy as np
+from sklearn.linear_model import LogisticRegression
+
 import quapy as qp
+import quapy.functional as F
 from quapy.data import LabelledCollection
 from quapy.functional import strprev
-from sklearn.linear_model import LogisticRegression
-import numpy as np
 from quapy.method.aggregative import PACC
-import quapy.functional as F
+from quapy.tests._synthetic import make_dataset
 
 
 class TestReplicability(unittest.TestCase):
 
-    def test_prediction_replicability(self):
+    @classmethod
+    def setUpClass(cls):
+        cls.binary_dataset = make_dataset(
+            n_train=180, n_test=80, n_classes=2, n_features=10, random_state=21, name='rep-binary'
+        )
+        cls.multiclass_dataset = make_dataset(
+            n_train=180, n_test=80, n_classes=3, n_features=12, random_state=23, name='rep-multiclass'
+        )
 
-        dataset = qp.datasets.fetch_UCIBinaryDataset('yeast')
-        train, test = dataset.train_test
+    def test_prediction_replicability(self):
+        train, test = self.binary_dataset.train_test
 
         with qp.util.temp_seed(0):
             lr = LogisticRegression(random_state=0, max_iter=10000)
@@ -29,7 +39,6 @@ class TestReplicability(unittest.TestCase):
 
         self.assertEqual(str_prev1, str_prev2)
 
-
     def test_samping_replicability(self):
 
         def equal_collections(c1, c2, value=True):
@@ -60,53 +69,33 @@ class TestReplicability(unittest.TestCase):
             sample2 = data.sampling(50, *[0.7, 0.3])
         equal_collections(sample1, sample2, True)
 
-        sample1 = data.sampling(50, *[0.7, 0.3], random_state=0)
-        sample2 = data.sampling(50, *[0.7, 0.3], random_state=0)
-        equal_collections(sample1, sample2, True)
-
         sample1_tr, sample1_te = data.split_stratified(train_prop=0.7, random_state=0)
         sample2_tr, sample2_te = data.split_stratified(train_prop=0.7, random_state=0)
         equal_collections(sample1_tr, sample2_tr, True)
         equal_collections(sample1_te, sample2_te, True)
 
-        with qp.util.temp_seed(0):
-            sample1_tr, sample1_te = data.split_stratified(train_prop=0.7)
-        with qp.util.temp_seed(0):
-            sample2_tr, sample2_te = data.split_stratified(train_prop=0.7)
-        equal_collections(sample1_tr, sample2_tr, True)
-        equal_collections(sample1_te, sample2_te, True)
-
-
     def test_parallel_replicability(self):
-
-        train, test = qp.datasets.fetch_UCIMulticlassDataset('dry-bean').reduce().train_test
-
-        test = test.sampling(500, *[0.1, 0.0, 0.1, 0.1, 0.2, 0.5, 0.0])
+        train, test = self.multiclass_dataset.train_test
+        test = test.sampling(60, *[0.2, 0.3, 0.5], random_state=4)
 
         with qp.util.temp_seed(10):
-            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
+            pacc = PACC(LogisticRegression(max_iter=5000), val_split=.5, n_jobs=2)
             pacc.fit(*train.Xy)
             prev1 = F.strprev(pacc.predict(test.instances))
 
         with qp.util.temp_seed(0):
-            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
+            pacc = PACC(LogisticRegression(max_iter=5000), val_split=.5, n_jobs=2)
             pacc.fit(*train.Xy)
             prev2 = F.strprev(pacc.predict(test.instances))
 
         with qp.util.temp_seed(0):
-            pacc = PACC(LogisticRegression(), val_split=.5, n_jobs=2)
+            pacc = PACC(LogisticRegression(max_iter=5000), val_split=.5, n_jobs=2)
             pacc.fit(*train.Xy)
             prev3 = F.strprev(pacc.predict(test.instances))
 
-        print(prev1)
-        print(prev2)
-        print(prev3)
-
         self.assertNotEqual(prev1, prev2)
         self.assertEqual(prev2, prev3)
 
 
-
-
 if __name__ == '__main__':
     unittest.main()