Add calibrated global min-PTE threshold for many null tests

src/sp_validation/statistics.py

@@ -3,10 +3,12 @@
 :Name: statistics.py
 
 :Description: Cosmology-independent statistical helpers (jackknife resampling,
-              chi2/PTE, covariance<->correlation, OneCovariance reshaping).
-              Extracted verbatim from the former basic.py.
+              chi2/PTE, calibrated min-PTE across many null tests,
+              covariance<->correlation, OneCovariance reshaping).
 """
 
+from dataclasses import dataclass
+
 import numpy as np
 from scipy import stats
 
@@ -126,3 +128,154 @@ def cov_from_one_covariance(cov_one_cov, gaussian=True):
         for j in range(n_bins):
             cov[i, j] = cov_one_cov[i * n_bins + j, index_value]
     return cov
+
+
+def _min_pte(ptes, two_sided):
+    """Minimum PTE across the last axis, optionally two-sided."""
+    ptes = np.asarray(ptes, dtype=float)
+    if not np.all((ptes >= 0.0) & (ptes <= 1.0)):
+        raise ValueError("PTEs must be finite and lie in [0, 1]")
+    if two_sided:
+        ptes = 2.0 * np.minimum(ptes, 1.0 - ptes)
+    return ptes.min(axis=-1)
+
+
+def wilson_interval(count, n, level=0.68):
+    """Two-sided Wilson score interval for a binomial fraction ``count / n``."""
+    z = stats.norm.ppf(0.5 + level / 2.0)
+    fraction = count / n
+    denominator = 1.0 + z**2 / n
+    centre = (fraction + z**2 / (2.0 * n)) / denominator
+    half_width = (
+        z * np.sqrt(fraction * (1.0 - fraction) / n + z**2 / (4.0 * n**2)) / denominator
+    )
+    return max(0.0, centre - half_width), min(1.0, centre + half_width)
+
+
+def effective_number_of_tests(threshold, alpha):
+    """Number of independent tests ``k`` with ``1 - (1 - threshold)^k = alpha``.
+
+    For ``k`` independent uniform PTEs, the ``alpha``-quantile of their minimum
+    is ``1 - (1 - alpha)^(1/k)``; inverting that for a calibrated threshold
+    gives the effective number of independent tests it corresponds to.
+    """
+    return np.log1p(-alpha) / np.log1p(-np.asarray(threshold, dtype=float))
+
+
+@dataclass(frozen=True)
+class MinPTECalibration:
+    """Global min-PTE null-test threshold calibrated on noise-only mocks.
+
+    Attributes
+    ----------
+    alpha : float
+        Global false-positive rate the threshold is calibrated to.
+    two_sided : bool
+        Whether each PTE ``p`` entered as ``2 min(p, 1 - p)``.
+    mock_min_pte : numpy.ndarray
+        Minimum PTE across statistics for each mock, shape ``(n_mocks,)``.
+    threshold : float
+        ``alpha``-quantile of ``mock_min_pte``: a data vector whose minimum
+        PTE falls below it fails the global null test at level ``alpha``.
+    threshold_interval : tuple of float
+        Distribution-free interval on ``threshold`` at confidence ``level``.
+    k_eff : float
+        Effective number of independent tests implied by ``threshold``.
+    k_eff_interval : tuple of float
+        ``k_eff`` mapped from ``threshold_interval``.
+    level : float
+        Confidence level of the intervals.
+    """
+
+    alpha: float
+    two_sided: bool
+    mock_min_pte: np.ndarray
+    threshold: float
+    threshold_interval: tuple
+    k_eff: float
+    k_eff_interval: tuple
+    level: float
+
+    @property
+    def n_mocks(self):
+        """Number of mocks the calibration rests on."""
+        return self.mock_min_pte.size
+
+    def global_pte(self, data_ptes):
+        """Global p-value of a data vector's PTEs against the mocks.
+
+        Parameters
+        ----------
+        data_ptes : array_like
+            The data's PTE for each statistic, ordered as the mock columns.
+
+        Returns
+        -------
+        tuple
+            ``(p, interval)``: the fraction of mocks whose minimum PTE is at
+            most the data's, and its Wilson interval at ``level``.
+        """
+        data_min = _min_pte(data_ptes, self.two_sided)
+        count = int(np.count_nonzero(self.mock_min_pte <= data_min))
+        return count / self.n_mocks, wilson_interval(count, self.n_mocks, self.level)
+
+
+def calibrate_min_pte(mock_ptes, alpha=0.05, two_sided=False, level=0.68):
+    """Calibrate a global threshold on the minimum PTE across statistics.
+
+    Testing many correlated statistics at a fixed per-test level inflates the
+    false-positive rate by an unknown amount. Taking the minimum PTE across
+    statistics for each noise-only mock and reading off its ``alpha``-quantile
+    gives a threshold with exactly that global rate, whatever the correlations
+    between statistics and whatever miscalibration of the individual PTEs.
+
+    Parameters
+    ----------
+    mock_ptes : array_like
+        PTEs of noise-only realisations, shape ``(n_mocks, n_stats)``; the
+        columns can be any statistics, redshift-bin pairs or scale cuts.
+    alpha : float, optional
+        Global false-positive rate; default is ``0.05``.
+    two_sided : bool, optional
+        If ``True``, test each PTE ``p`` as ``2 min(p, 1 - p)`` so that
+        anomalously small statistics also fail; default is ``False``.
+    level : float, optional
+        Confidence level of the reported intervals; default is ``0.68``.
+
+    Returns
+    -------
+    MinPTECalibration
+        The threshold, its interval, the effective number of independent
+        tests, and :meth:`MinPTECalibration.global_pte` for the data.
+
+    Notes
+    -----
+    The interval on the threshold uses order statistics: the number of mocks
+    below the true ``alpha``-quantile is ``Binomial(n_mocks, alpha)``, and the
+    ``(1 -/+ level) / 2`` quantiles of that count index the mocks bracketing it.
+    """
+    mock_ptes = np.asarray(mock_ptes, dtype=float)
+    if mock_ptes.ndim != 2:
+        raise ValueError("mock_ptes must have shape (n_mocks, n_stats)")
+    if not 0.0 < alpha < 1.0:
+        raise ValueError("alpha must lie in (0, 1)")
+
+    minima = _min_pte(mock_ptes, two_sided)
+    n = minima.size
+    ordered = np.sort(minima)
+    tail = (1.0 - level) / 2.0
+    lower = int(np.clip(stats.binom.ppf(tail, n, alpha), 1, n)) - 1
+    upper = int(np.clip(stats.binom.ppf(1.0 - tail, n, alpha) + 1, 1, n)) - 1
+    threshold = float(np.quantile(minima, alpha))
+    interval = (float(ordered[lower]), float(ordered[upper]))
+    k_eff = effective_number_of_tests([threshold, *interval], alpha)
+    return MinPTECalibration(
+        alpha=alpha,
+        two_sided=two_sided,
+        mock_min_pte=minima,
+        threshold=threshold,
+        threshold_interval=interval,
+        k_eff=float(k_eff[0]),
+        k_eff_interval=(float(k_eff[2]), float(k_eff[1])),
+        level=level,
+    )

src/sp_validation/tests/test_statistics.py

@@ -15,9 +15,11 @@
 from scipy import stats
 
 from sp_validation.statistics import (
+    calibrate_min_pte,
     chi2_and_pte,
     corr_from_cov,
     cov_from_one_covariance,
+    effective_number_of_tests,
     jackknif_weighted_average2,
 )
 
@@ -194,3 +196,60 @@ def test_cov_from_one_covariance_selects_gaussian_column():
         [[0.0, 1.0], [10.0, 11.0]],
         rtol=1e-12,
     )
+
+
+def test_calibrate_min_pte_independent_uniform_ptes():
+    """Independent uniform PTEs recover the Sidak threshold and k_eff = k."""
+    rng = np.random.default_rng(1)
+    alpha, k = 0.05, 6
+    cal = calibrate_min_pte(rng.uniform(size=(200_000, k)), alpha=alpha)
+    expected = 1.0 - (1.0 - alpha) ** (1.0 / k)
+    npt.assert_allclose(cal.threshold, expected, rtol=0.02)
+    npt.assert_allclose(cal.k_eff, k, rtol=0.02)
+    assert cal.threshold_interval[0] <= cal.threshold <= cal.threshold_interval[1]
+    assert cal.k_eff_interval[0] <= cal.k_eff <= cal.k_eff_interval[1]
+    npt.assert_allclose(effective_number_of_tests(expected, alpha), k)
+
+
+def test_calibrate_min_pte_perfectly_correlated_is_one_test():
+    """Identical columns collapse to one test: threshold alpha, k_eff 1."""
+    rng = np.random.default_rng(2)
+    column = rng.uniform(size=(100_000, 1))
+    cal = calibrate_min_pte(np.repeat(column, 8, axis=1), alpha=0.05)
+    npt.assert_allclose(cal.threshold, 0.05, rtol=0.03)
+    npt.assert_allclose(cal.k_eff, 1.0, rtol=0.03)
+
+
+def test_calibrate_min_pte_two_sided_independent():
+    """Two-sided PTEs 2 min(p, 1 - p) are uniform, so Sidak still holds."""
+    rng = np.random.default_rng(3)
+    cal = calibrate_min_pte(rng.uniform(size=(200_000, 4)), alpha=0.05, two_sided=True)
+    npt.assert_allclose(cal.k_eff, 4.0, rtol=0.03)
+
+
+def test_min_pte_global_pte():
+    """Global p-value is the mock fraction with min PTE <= the data's."""
+    mock_ptes = np.array([[0.1, 0.9], [0.5, 0.2], [0.3, 0.7], [0.8, 0.6]])
+    cal = calibrate_min_pte(mock_ptes, alpha=0.25)
+    # Mock minima: 0.1, 0.2, 0.3, 0.6.
+    p, (lo, hi) = cal.global_pte([0.9, 0.2])
+    assert p == 0.5
+    assert lo < 0.5 < hi
+    assert cal.global_pte([0.05, 0.5])[0] == 0.0
+    assert cal.global_pte([0.99, 0.95])[0] == 1.0
+    # Two-sided: a suspiciously good PTE of 0.99 counts like 0.02.
+    two = calibrate_min_pte(mock_ptes, alpha=0.25, two_sided=True)
+    assert two.global_pte([0.99, 0.5])[0] == 0.0
+
+
+def test_global_pte_is_calibrated_under_the_null():
+    """For null data the global p-value is uniform: P(p <= alpha) ~ alpha."""
+    rng = np.random.default_rng(4)
+    mean = np.zeros(5)
+    cov = 0.6 * np.ones((5, 5)) + 0.4 * np.eye(5)
+    to_pte = lambda z: stats.norm.sf(z)  # noqa: E731
+    cal = calibrate_min_pte(to_pte(rng.multivariate_normal(mean, cov, 4000)))
+    data = to_pte(rng.multivariate_normal(mean, cov, 4000))
+    p = np.array([cal.global_pte(row)[0] for row in data])
+    npt.assert_allclose(np.mean(p <= 0.05), 0.05, atol=0.012)
+    assert 1.0 < cal.k_eff < 5.0