diff --git a/docs/api/calib/pyhealth.calib.predictionset.rst b/docs/api/calib/pyhealth.calib.predictionset.rst index 740b1c87c..9f00edebb 100644 --- a/docs/api/calib/pyhealth.calib.predictionset.rst +++ b/docs/api/calib/pyhealth.calib.predictionset.rst @@ -8,10 +8,14 @@ techniques for uncertainty quantification. ``BaseConformal``, ``LABEL``, ``ClusterLabel``, ``CovariateLabel``, and ``NeighborhoodLabel`` all accept a ``score_type`` argument selecting the nonconformity/conformity score used for calibration and set construction: -either ``"threshold"`` (the default, unchanged from prior releases) or -``"aps"`` (Adaptive Prediction Sets, Romano, Sesia, and Candes 2020), which -adapts the prediction set size to the model's per-input confidence. See -:mod:`pyhealth.calib.predictionset.scores` for the exact score formulas. +``"threshold"`` (LAC, Sadinle, Lei, and Wasserman 2019 -- the default for +every one of these classes) or ``"aps"`` (Adaptive Prediction Sets, Romano, +Sesia, and Candes 2020, which adapts the prediction set size to the +model's per-input confidence). ``BaseConformal`` additionally supports +``"margin"`` (Papadopoulos, Vovk, and Gammerman 2007's own score for +neural-network classifiers). Coverage validity does not depend on this +choice -- see :mod:`pyhealth.calib.predictionset.scores` for the exact +score formulas. ``SCRIB`` and ``FavMac`` are not included since their calibration procedures aren't a score-then-quantile pattern. @@ -30,8 +34,26 @@ Available Methods pyhealth.calib.predictionset.ClusterLabel pyhealth.calib.predictionset.NeighborhoodLabel -BaseConformal (Standard Split Conformal Prediction) ----------------------------------------------------- +BaseConformal (Standard "Naive" Split Conformal Prediction) +-------------------------------------------------------------- + +.. note:: + + ``BaseConformal`` and ``LABEL`` both build on the same generic + split-conformal calibration machinery (Vovk, Gammerman, and Shafer + 2005): sort calibration nonconformity scores and take the + ``ceil((1-alpha)*(n+1))``-th one as the threshold. By default they use + the identical score too -- ``"threshold"``, i.e. ``1 - p(true class)``, + from Sadinle, Lei, and Wasserman (2019) -- for consistency with every + other class in this module. ``BaseConformal`` additionally supports + ``score_type="margin"``, i.e. ``max_{j!=k} p(j) - p(k)``, the + nonconformity measure from Papadopoulos, Vovk, and Gammerman, + "Conformal prediction with neural networks," 19th IEEE ICTAI 2007, + vol. 2, pp. 388-395 (verified directly against Section 4.2 of + Papadopoulos's 2008 InTech chapter restating that paper, which calls + this "the natural nonconformity measure" for neural-network + classifiers). Coverage validity does not depend on this choice; it + only changes prediction-set composition and average set size. .. autoclass:: pyhealth.calib.predictionset.BaseConformal :members: diff --git a/examples/conformal_eeg/tuev_conventional_conformal.py b/examples/conformal_eeg/tuev_conventional_conformal.py index e5542c7d4..11f842115 100644 --- a/examples/conformal_eeg/tuev_conventional_conformal.py +++ b/examples/conformal_eeg/tuev_conventional_conformal.py @@ -17,6 +17,16 @@ python examples/conformal_eeg/tuev_conventional_conformal.py \\ --root downloads/tuev/v2.0.1/edf --model tfm --n-seeds 5 --seed 42 --alpha 0.1 \\ --log-file tuev_conventional_tfm_alpha0.1_5seeds.log + +Note: this example specifically demonstrates LABEL (Sadinle, Lei, and +Wasserman 2019's least-ambiguous-set score, "1 - p(true class)"), which is +also pyhealth.calib.predictionset.BaseConformal's default score_type +(both classes default to "threshold" for consistency). BaseConformal +additionally supports score_type="margin", the nonconformity measure from +Papadopoulos, Vovk, and Gammerman's "Conformal prediction with neural +networks" (ICTAI 2007) -- pass BaseConformal(model, alpha, +score_type="margin") if you specifically want that paper's own method. +Coverage validity does not depend on which score_type is used. """ from __future__ import annotations diff --git a/pyhealth/calib/predictionset/base_conformal/__init__.py b/pyhealth/calib/predictionset/base_conformal/__init__.py index 9dde35db5..bb5d69da3 100644 --- a/pyhealth/calib/predictionset/base_conformal/__init__.py +++ b/pyhealth/calib/predictionset/base_conformal/__init__.py @@ -1,8 +1,19 @@ """ Base Conformal Prediction (Split Conformal) -Standard split conformal prediction for multiclass classification without -covariate shift correction. +Standard "naive" split conformal prediction for multiclass classification +without covariate shift correction. Like every other score-then-quantile +class in this package (LABEL, ClusterLabel, CovariateLabel, +NeighborhoodLabel), the default nonconformity score is "threshold" +(Sadinle, Lei, and Wasserman 2019's least-ambiguous-set score), for +consistency across the module. BaseConformal additionally supports +score_type="margin", the nonconformity measure Papadopoulos, Vovk, and +Gammerman defined specifically for neural-network classifiers -- pass this +if you want that paper's own method rather than LABEL's score. Coverage +validity does not depend on which score_type is chosen (Vovk, Gammerman, +and Shafer 2005): all of them share the same split-conformal calibration +machinery and guarantee, and differ only in prediction-set composition +and average set size. This method constructs prediction sets with coverage guarantees by calibrating score thresholds on a held-out calibration set. @@ -11,12 +22,15 @@ Vovk, Vladimir, Alexander Gammerman, and Glenn Shafer. "Algorithmic learning in a random world." Springer, 2005. - Papadopoulos, Harris, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. - "Inductive confidence machines for regression." ECML 2002. - Sadinle, Mauricio, Jing Lei, and Larry Wasserman. "Least ambiguous set-valued classifiers with bounded error levels." Journal of the - American Statistical Association (2019). [score_type="threshold"] + American Statistical Association (2019). [score_type="threshold", + the default] + + Papadopoulos, Harris, Vladimir Vovk, and Alexander Gammerman. + "Conformal prediction with neural networks." 19th IEEE International + Conference on Tools with Artificial Intelligence (ICTAI 2007), vol. 2, + pp. 388-395. IEEE, 2007. [score_type="margin"] Romano, Yaniv, Matteo Sesia, and Emmanuel Candes. "Classification with valid and adaptive coverage." NeurIPS 2020. [score_type="aps"] @@ -95,26 +109,45 @@ def _query_weighted_quantile( class BaseConformal(SetPredictor): - """Base Conformal Prediction for multiclass classification. + """Base ("naive") Conformal Prediction for multiclass classification. This implements standard split conformal prediction, which constructs prediction sets with distribution-free coverage guarantees. The method calibrates thresholds on a calibration set and uses them to construct prediction sets on test data. + By default this uses the same "threshold" nonconformity score as every + other class in this package (LABEL, ClusterLabel, CovariateLabel, + NeighborhoodLabel): Sadinle, Lei, and Wasserman's (2019) ``1 - p(k)`` + score. BaseConformal additionally supports score_type="margin", the + nonconformity measure Papadopoulos, Vovk, and Gammerman (2007) defined + for neural-network classifiers: ``alpha(x, k) = max_{j != k} p(j) - + p(k)``, i.e. how much the strongest competing class beats k. Coverage + validity does not depend on this choice (Vovk, Gammerman, and Shafer + 2005): any nonconformity score, under the same rank-based calibration + procedure, gives the same coverage guarantee -- score_type only changes + prediction-set composition and average set size, not validity. + The method guarantees that: - For marginal coverage (alpha is float): P(Y not in C(X)) <= alpha - For class-conditional coverage (alpha is array): P(Y not in C(X) | Y=k) <= alpha[k] - where C(X) denotes the prediction set for input X. + where C(X) denotes the prediction set for input X. This holds regardless + of score_type, since validity of split conformal prediction does not + depend on which nonconformity measure is used (only the usefulness / + average set size does). Papers: Vovk, Vladimir, Alexander Gammerman, and Glenn Shafer. "Algorithmic learning in a random world." Springer, 2005. - Lei, Jing, Max G'Sell, Alessandro Rinaldo, Ryan J. Tibshirani, - and Larry Wasserman. "Distribution-free predictive inference for - regression." Journal of the American Statistical Association (2018). + Sadinle, Mauricio, Jing Lei, and Larry Wasserman. "Least ambiguous + set-valued classifiers with bounded error levels." Journal of the + American Statistical Association (2019). + + Papadopoulos, Harris, Vladimir Vovk, and Alexander Gammerman. + "Conformal prediction with neural networks." 19th IEEE ICTAI 2007, + vol. 2, pp. 388-395. Args: model: A trained base model that outputs predicted probabilities @@ -122,8 +155,13 @@ class BaseConformal(SetPredictor): - float: marginal coverage P(Y not in C(X)) <= alpha - array: class-conditional P(Y not in C(X) | Y=k) <= alpha[k] score_type: Type of nonconformity score to use: - - "threshold" (default): NC score = 1 - p(true class), the score - from Sadinle, Lei, and Wasserman (2019) ("LABEL"). + - "threshold" (default): NC score = 1 - p(true class), the + score from Sadinle, Lei, and Wasserman (2019) ("LABEL") -- + matches :class:`~pyhealth.calib.predictionset.LABEL`'s + behavior, and every other class in this package's default. + - "margin": NC score = max_{j!=k} p(j) - p(k), the score from + Papadopoulos, Vovk, and Gammerman (2007) ("naive" split + conformal prediction for neural networks). - "aps": Adaptive Prediction Sets (Romano, Sesia, and Candes 2020). NC score for class k is the cumulative sum of predicted probabilities for classes ranked above k, plus a randomized @@ -134,7 +172,8 @@ class BaseConformal(SetPredictor): individual input. See :mod:`pyhealth.calib.predictionset.scores` for the exact formula. random_state: Optional int seed for the RNG used by score_type="aps" - (the U ~ Uniform(0,1) draws). Ignored for score_type="threshold". + (the U ~ Uniform(0,1) draws). Ignored for score_type="threshold" + and score_type="margin". debug: Whether to use debug mode (processes fewer samples) Examples: @@ -181,6 +220,12 @@ class BaseConformal(SetPredictor): >>> conformal_model_aps = BaseConformal( ... model, alpha=0.1, score_type="aps", random_state=0) >>> conformal_model_aps.calibrate(cal_dataset=val_data) + >>> + >>> # Use Papadopoulos, Vovk, and Gammerman (2007)'s own margin + >>> # score instead of the default threshold (LABEL) score + >>> conformal_model_margin = BaseConformal( + ... model, alpha=0.1, score_type="margin") + >>> conformal_model_margin.calibrate(cal_dataset=val_data) """ def __init__( diff --git a/pyhealth/calib/predictionset/scores.py b/pyhealth/calib/predictionset/scores.py index ae4172d1d..dd5e6b66f 100644 --- a/pyhealth/calib/predictionset/scores.py +++ b/pyhealth/calib/predictionset/scores.py @@ -2,10 +2,11 @@ This module separates the *score* used by a conformal-prediction-set method from the *calibration/thresholding procedure* it's plugged into. These are -two independent axes: the choice of score ("threshold"/LAC vs "aps") does -not depend on how the resulting scores get turned into a threshold (marginal -quantile, per-class quantile, per-cluster quantile, weighted quantile for -covariate shift, or localized weighted quantile for neighborhood methods). +two independent axes: the choice of score ("threshold"/LAC vs "aps" vs +"margin") does not depend on how the resulting scores get turned into a +threshold (marginal quantile, per-class quantile, per-cluster quantile, +weighted quantile for covariate shift, or localized weighted quantile for +neighborhood methods). Supported score types: @@ -31,9 +32,26 @@ to the model's confidence for each individual input, which the "threshold" score does not. -Both scores are computed here in *nonconformity* convention (higher = less -conforming, i.e. 1 minus a probability-like quantity) since that's the -convention BaseConformal/LABEL/ClusterLabel use internally. A *conformity* + - "margin" (Papadopoulos, Vovk, and Gammerman, "Conformal Prediction + with Neural Networks," 19th IEEE ICTAI 2007, vol. 2, pp. 388-395; + restated in Section 4.2 of Papadopoulos, "Inductive Conformal + Prediction: Theory and Application to Neural Networks," Tools in + Artificial Intelligence, InTech, 2008, as the "natural nonconformity + measure" for neural-network classifiers): the score for class k is + how much the best *other* class beats k:: + + alpha(x, k) = max_{j != k} pi(x, j) - pi(x, k) + + Unlike "threshold" (which only looks at k's own probability), this + measure the margin between k and its strongest competitor, so it is + more nonconforming for an example whose true class is only narrowly + ahead of a rival than for one that's ahead by a landslide, even at + the same raw probability for k. + +Both "threshold" and "aps" scores are computed here in *nonconformity* +convention (higher = less conforming, i.e. 1 minus a probability-like +quantity) since that's the convention BaseConformal/LABEL/ClusterLabel use +internally; "margin" is nonconformity-signed by construction. A *conformity* (higher = more conforming) variant is also provided for CovariateLabel/ NeighborhoodLabel, which use the opposite sign convention internally; it is simply `1 - nonconformity`, preserving the same ranking of examples either @@ -50,7 +68,7 @@ "true_class_nc_scores", ] -SUPPORTED_SCORE_TYPES = ("threshold", "aps") +SUPPORTED_SCORE_TYPES = ("threshold", "aps", "margin") def _validate_score_type(score_type: str) -> None: @@ -112,6 +130,38 @@ def _aps_all_class_nc_scores( return scores +def _margin_all_class_nc_scores(y_prob: np.ndarray) -> np.ndarray: + """Computes the margin nonconformity score for every class, every row. + + alpha(x, k) = max_{j != k} pi(x, j) - pi(x, k) + + (Papadopoulos, Vovk, and Gammerman 2007, Eq. for the "natural + nonconformity measure"; see module docstring.) Vectorized via each + row's top-2 probabilities, with ties at the row max handled so that a + class tied for the top spot still has another class achieving that same + max value as its "best other class." + + Args: + y_prob: Predicted probabilities, shape (N, K). + + Returns: + Nonconformity scores of shape (N, K); higher means less conforming. + """ + n, k = y_prob.shape + order = np.argsort(-y_prob, axis=1) + sorted_probs = np.take_along_axis(y_prob, order, axis=1) + top1 = sorted_probs[:, 0] + top2 = sorted_probs[:, 1] if k >= 2 else np.full(n, -np.inf) + is_top = y_prob == top1[:, None] + num_top = is_top.sum(axis=1, keepdims=True) + max_others = np.where( + is_top, + np.where(num_top > 1, top1[:, None], top2[:, None]), + top1[:, None], + ) + return max_others - y_prob + + def all_class_nc_scores( y_prob: np.ndarray, score_type: str = "threshold", @@ -122,12 +172,13 @@ def all_class_nc_scores( Args: y_prob: Predicted probabilities, shape (N, K). - score_type: "threshold" (Sadinle, Lei, and Wasserman 2019) or "aps" - (Romano, Sesia, and Candes 2020). Default "threshold". + score_type: "threshold" (Sadinle, Lei, and Wasserman 2019), "aps" + (Romano, Sesia, and Candes 2020), or "margin" (Papadopoulos, + Vovk, and Gammerman 2007). Default "threshold". rng: Random generator, required (and only used) if score_type="aps" and randomize=True. randomize: Whether to use the randomized ("exact coverage") variant - of APS. Ignored for score_type="threshold". + of APS. Ignored for score_type="threshold" and "margin". Returns: Nonconformity scores of shape (N, K). @@ -144,10 +195,15 @@ def all_class_nc_scores( >>> np.round(scores, 2) array([[0.42, 0.82, 0.96], [0.72, 0.36, 0.95]]) + >>> all_class_nc_scores(y_prob, score_type="margin") + array([[-0.5, 0.5, 0.6], + [ 0.2, -0.2, 0.3]]) """ _validate_score_type(score_type) if score_type == "threshold": return 1.0 - y_prob + if score_type == "margin": + return _margin_all_class_nc_scores(y_prob) # score_type == "aps" if rng is None: rng = np.random.default_rng() diff --git a/tests/core/test_base_conformal.py b/tests/core/test_base_conformal.py new file mode 100644 index 000000000..95972d23d --- /dev/null +++ b/tests/core/test_base_conformal.py @@ -0,0 +1,269 @@ +"""Tests for BaseConformal, standard split conformal prediction for +multiclass classification (Vovk, Gammerman, and Shafer 2005). + +BaseConformal's default score_type ("threshold") is Sadinle, Lei, and +Wasserman's (2019) LAC score, matching every other score-then-quantile +class in this package (LABEL, ClusterLabel, CovariateLabel, +NeighborhoodLabel) for consistency. BaseConformal additionally supports +score_type="margin", the nonconformity measure Papadopoulos, Vovk, and +Gammerman defined for neural-network classifiers ("Conformal prediction +with neural networks," 19th IEEE ICTAI 2007, vol. 2, pp. 388-395) for +anyone who specifically wants that paper's own method. These tests verify +both score types work correctly, that coverage validity holds for both +(as it must for any nonconformity score), and that "margin" is a real, +distinct option rather than a no-op alias for "threshold". +""" + +import unittest + +import numpy as np +import torch + +from pyhealth.calib.predictionset.base_conformal import BaseConformal +from pyhealth.calib.predictionset.label import LABEL +from pyhealth.datasets import create_sample_dataset, get_dataloader +from pyhealth.models import MLP + + +class TestBaseConformal(unittest.TestCase): + """Test cases for the BaseConformal prediction set constructor.""" + + def setUp(self): + np.random.seed(42) + torch.manual_seed(42) + + # 3-class multiclass task, split into a 6-sample train set (indices + # 0-5) and a 12-sample calibration set (indices 6-17) so quantile + # thresholds are well-defined even at small alpha. + self.samples = [ + { + "patient_id": f"patient-{i}", + "visit_id": f"visit-{i}", + "conditions": [f"cond-{i}", f"cond-{i+1}", f"cond-{i+2}"], + "procedures": [float(i), float(i + 1), float(i + 2), float(i + 3)], + "label": i % 3, + } + for i in range(18) + ] + + self.input_schema = {"conditions": "sequence", "procedures": "tensor"} + self.output_schema = {"label": "multiclass"} + + self.dataset = create_sample_dataset( + samples=self.samples, + input_schema=self.input_schema, + output_schema=self.output_schema, + dataset_name="test", + ) + + self.model = MLP( + dataset=self.dataset, + feature_keys=["conditions", "procedures"], + label_key="label", + mode="multiclass", + ) + self.model.eval() + + self.train_indices = list(range(6)) + self.cal_indices = list(range(6, 18)) + self.cal_dataset = self.dataset.subset(self.cal_indices) + + # -- initialization -------------------------------------------------- + + def test_default_score_type_is_threshold(self): + """BaseConformal's default must be "threshold" (Sadinle, Lei, and + Wasserman 2019's LAC score), matching every other score-then- + quantile class in this package (LABEL, ClusterLabel, + CovariateLabel, NeighborhoodLabel) for consistency.""" + base_model = BaseConformal(model=self.model, alpha=0.1) + self.assertEqual(base_model.score_type, "threshold") + + def test_initialization_with_array_alpha(self): + alpha_per_class = [0.1, 0.15, 0.2] + base_model = BaseConformal(model=self.model, alpha=alpha_per_class) + self.assertIsInstance(base_model.alpha, np.ndarray) + np.testing.assert_array_equal(base_model.alpha, alpha_per_class) + + def test_initialization_non_multiclass_raises_error(self): + binary_samples = [ + { + "patient_id": "patient-0", + "visit_id": "visit-0", + "conditions": ["cond-1"], + "procedures": [1.0], + "label": 0, + }, + { + "patient_id": "patient-1", + "visit_id": "visit-1", + "conditions": ["cond-2"], + "procedures": [2.0], + "label": 1, + }, + ] + binary_dataset = create_sample_dataset( + samples=binary_samples, + input_schema={"conditions": "sequence", "procedures": "tensor"}, + output_schema={"label": "binary"}, + dataset_name="test", + ) + binary_model = MLP( + dataset=binary_dataset, + feature_keys=["conditions"], + label_key="label", + mode="binary", + ) + with self.assertRaises(NotImplementedError): + BaseConformal(model=binary_model, alpha=0.1) + + def test_invalid_score_type_raises(self): + with self.assertRaises(ValueError): + BaseConformal(model=self.model, alpha=0.1, score_type="not_a_score") + + # -- calibration / forward, per score_type ---------------------------- + + def test_calibrate_and_forward_marginal_default_score(self): + base_model = BaseConformal(model=self.model, alpha=0.3) + base_model.calibrate(cal_dataset=self.cal_dataset) + self.assertIsNotNone(base_model.t) + # Marginal coverage -> a single scalar threshold. + self.assertEqual(base_model.t.numel(), 1) + + test_loader = get_dataloader(self.dataset, batch_size=2, shuffle=False) + data_batch = next(iter(test_loader)) + with torch.no_grad(): + output = base_model(**data_batch) + self.assertIn("y_predset", output) + self.assertEqual(output["y_predset"].dtype, torch.bool) + self.assertEqual(output["y_predset"].shape, output["y_prob"].shape) + + def test_calibrate_class_conditional(self): + alpha_per_class = [0.3, 0.35, 0.3] + base_model = BaseConformal(model=self.model, alpha=alpha_per_class) + base_model.calibrate(cal_dataset=self.cal_dataset) + self.assertEqual(base_model.t.numel(), 3) + + def test_score_type_threshold_runs_end_to_end(self): + base_model = BaseConformal(model=self.model, alpha=0.3, score_type="threshold") + base_model.calibrate(cal_dataset=self.cal_dataset) + test_loader = get_dataloader(self.dataset, batch_size=2, shuffle=False) + with torch.no_grad(): + for data_batch in test_loader: + output = base_model(**data_batch) + self.assertEqual(output["y_predset"].dtype, torch.bool) + set_sizes = output["y_predset"].sum(dim=1) + self.assertTrue(torch.all(set_sizes > 0)) + + def test_score_type_aps_runs_end_to_end(self): + base_model = BaseConformal( + model=self.model, alpha=0.3, score_type="aps", random_state=42 + ) + base_model.calibrate(cal_dataset=self.cal_dataset) + test_loader = get_dataloader(self.dataset, batch_size=2, shuffle=False) + with torch.no_grad(): + for data_batch in test_loader: + output = base_model(**data_batch) + self.assertEqual(output["y_predset"].dtype, torch.bool) + set_sizes = output["y_predset"].sum(dim=1) + self.assertTrue(torch.all(set_sizes > 0)) + + def test_forward_before_calibration_raises_error(self): + base_model = BaseConformal(model=self.model, alpha=0.2) + test_loader = get_dataloader(self.dataset, batch_size=1, shuffle=False) + data_batch = next(iter(test_loader)) + with self.assertRaises(RuntimeError): + with torch.no_grad(): + base_model(**data_batch) + + # -- the core claim: BaseConformal now genuinely differs from LABEL --- + + def test_default_matches_label_default(self): + """BaseConformal's default (score_type="threshold") should exactly + reproduce LABEL's thresholds and prediction sets on the same data, + since both then use the identical Sadinle et al. 2019 score and + calibration procedure -- this is intentional consistency across + the package's CP classes, not a bug.""" + alpha = 0.3 + base_model = BaseConformal(model=self.model, alpha=alpha) + base_model.calibrate(cal_dataset=self.cal_dataset) + + label_model = LABEL(model=self.model, alpha=alpha) + label_model.calibrate(cal_dataset=self.cal_dataset) + + self.assertAlmostEqual( + float(base_model.t.item()), float(label_model.t.item()), places=6 + ) + + test_loader = get_dataloader(self.dataset, batch_size=18, shuffle=False) + data_batch = next(iter(test_loader)) + with torch.no_grad(): + base_out = base_model(**data_batch) + label_out = label_model(**data_batch) + self.assertTrue(torch.equal(base_out["y_predset"], label_out["y_predset"])) + + def test_margin_score_type_diverges_from_default(self): + """score_type="margin" (Papadopoulos, Vovk, and Gammerman 2007's + own nonconformity measure) must produce a genuinely different + threshold and prediction sets than the default "threshold" score + on the same calibration data -- proving "margin" is a real, + distinct, working option and not a no-op alias.""" + alpha = 0.3 + default_model = BaseConformal(model=self.model, alpha=alpha) + default_model.calibrate(cal_dataset=self.cal_dataset) + + margin_model = BaseConformal(model=self.model, alpha=alpha, score_type="margin") + margin_model.calibrate(cal_dataset=self.cal_dataset) + + self.assertNotAlmostEqual( + float(default_model.t.item()), + float(margin_model.t.item()), + places=6, + msg="score_type=\"margin\" produced the same threshold as the " + "default \"threshold\" score -- expected a genuinely different " + "nonconformity measure.", + ) + + test_loader = get_dataloader(self.dataset, batch_size=18, shuffle=False) + data_batch = next(iter(test_loader)) + with torch.no_grad(): + default_out = default_model(**data_batch) + margin_out = margin_model(**data_batch) + + self.assertFalse( + torch.equal(default_out["y_predset"], margin_out["y_predset"]), + "score_type=\"margin\" produced identical prediction sets to " + "the default on every example.", + ) + + # -- coverage validity (holds regardless of score_type) --------------- + + def test_margin_score_achieves_approximate_marginal_coverage(self): + """Monte Carlo check: with a larger synthetic exchangeable dataset, + the "margin" score (available via score_type="margin") should + achieve close to the target 1-alpha marginal coverage, per split + conformal's validity guarantee (Vovk, Gammerman, and Shafer 2005), + which holds for any nonconformity measure.""" + rng = np.random.default_rng(123) + n, k = 3000, 4 + alpha = 0.1 + logits = rng.normal(size=(n, k)) * 2 + y_prob = np.exp(logits) / np.exp(logits).sum(1, keepdims=True) + y_true = np.array([rng.choice(k, p=y_prob[i]) for i in range(n)]) + + from pyhealth.calib.predictionset.base_conformal import _query_quantile + from pyhealth.calib.predictionset.scores import ( + all_class_nc_scores, + true_class_nc_scores, + ) + + cal, test = slice(0, n // 2), slice(n // 2, n) + nc_cal = true_class_nc_scores(y_prob[cal], y_true[cal], score_type="margin") + t = _query_quantile(nc_cal, alpha) + nc_test = all_class_nc_scores(y_prob[test], score_type="margin") + predset = nc_test <= t + coverage = predset[np.arange(n - n // 2), y_true[test]].mean() + self.assertGreaterEqual(coverage, 1 - alpha - 0.03) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/core/test_scores.py b/tests/core/test_scores.py index e6af15c69..c5c4dc359 100644 --- a/tests/core/test_scores.py +++ b/tests/core/test_scores.py @@ -1,6 +1,7 @@ """Tests for pyhealth.calib.predictionset.scores: the shared score module -implementing both the "threshold" (LAC) and "aps" (Adaptive Prediction -Sets, Romano/Sesia/Candes 2020) nonconformity/conformity scores. +implementing the "threshold" (LAC), "aps" (Adaptive Prediction Sets, +Romano/Sesia/Candes 2020), and "margin" (Papadopoulos/Vovk/Gammerman 2007) +nonconformity/conformity scores. """ import unittest @@ -112,6 +113,78 @@ def test_true_class_score_matches_all_class_indexing(self): np.testing.assert_allclose(true_scores, all_scores[np.arange(2), y_true]) +class TestScoresMargin(unittest.TestCase): + """"margin" is Papadopoulos, Vovk, and Gammerman (2007)'s nonconformity + measure for neural-network classifiers: alpha(x, k) = max_{j!=k} p(j) + - p(k). Verified against the ICTAI 2007 paper's Section 4.2 formula + (also restated verbatim in Papadopoulos's 2008 InTech chapter as "the + natural nonconformity measure").""" + + def test_hand_computed_no_ties(self): + y_prob = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]]) + # Row 0: top1=0.7 (class 0), top2=0.2. + # class0 (is top): 0.2 - 0.7 = -0.5 + # class1: 0.7 - 0.2 = 0.5 + # class2: 0.7 - 0.1 = 0.6 + # Row 1: top1=0.5 (class 1), top2=0.3. + # class0: 0.5 - 0.3 = 0.2 + # class1 (is top): 0.3 - 0.5 = -0.2 + # class2: 0.5 - 0.2 = 0.3 + scores = all_class_nc_scores(y_prob, score_type="margin") + np.testing.assert_allclose(scores, [[-0.5, 0.5, 0.6], [0.2, -0.2, 0.3]]) + + def test_true_class_matches_all_class_indexing(self): + y_prob = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]]) + y_true = np.array([0, 1]) + np.testing.assert_allclose( + true_class_nc_scores(y_prob, y_true, score_type="margin"), + [-0.5, -0.2], + ) + + def test_conformity_is_one_minus_nonconformity(self): + y_prob = np.array([[0.7, 0.2, 0.1], [0.3, 0.5, 0.2]]) + nc = all_class_nc_scores(y_prob, score_type="margin") + conf = all_class_conformity_scores(y_prob, score_type="margin") + np.testing.assert_allclose(conf, 1.0 - nc) + + def test_tie_at_max_uses_max_as_max_others_for_both_tied_classes(self): + """When two classes tie for the top probability, each tied class's + "best other class" is still the max value (the other tied class), + not the third-place value -- a naive "drop my own rank-1 slot" + implementation would get this wrong.""" + y_prob = np.array([[0.45, 0.45, 0.1]]) + scores = all_class_nc_scores(y_prob, score_type="margin")[0] + # Both tied classes: max_{j!=k} p_j = 0.45 (the other tied class). + np.testing.assert_allclose(scores[0], 0.45 - 0.45) + np.testing.assert_allclose(scores[1], 0.45 - 0.45) + # Untied third class: max_{j!=k} p_j = 0.45 (either tied class). + np.testing.assert_allclose(scores[2], 0.45 - 0.1) + + def test_true_class_with_highest_probability_has_lowest_score(self): + """The true class being far ahead of its best competitor should + yield a strongly negative (very conforming) margin score.""" + y_prob = np.array([[0.9, 0.06, 0.04]]) + y_true = np.array([0]) + score = true_class_nc_scores(y_prob, y_true, score_type="margin")[0] + self.assertLess(score, 0) + np.testing.assert_allclose(score, 0.06 - 0.9) + + def test_binary_classification_two_classes(self): + """With K=2, each class's only "other" is the remaining class.""" + y_prob = np.array([[0.8, 0.2], [0.35, 0.65]]) + scores = all_class_nc_scores(y_prob, score_type="margin") + np.testing.assert_allclose( + scores, [[0.2 - 0.8, 0.8 - 0.2], [0.65 - 0.35, 0.35 - 0.65]] + ) + + def test_scores_sum_to_zero_within_row_for_binary(self): + """For K=2 specifically, alpha(x,0) = -alpha(x,1) by construction.""" + rng = np.random.default_rng(8) + y_prob = rng.dirichlet([1, 1], size=20) + scores = all_class_nc_scores(y_prob, score_type="margin") + np.testing.assert_allclose(scores[:, 0], -scores[:, 1], atol=1e-12) + + class TestScoresCoverage(unittest.TestCase): """The core statistical property: both score types must achieve approximately the target marginal coverage under split conformal