Maxtest

src/pqm/__init__.py

@@ -1,5 +1,6 @@
 from .pqm import pqm_pvalue, pqm_chi2
 from .test_gaussian import test
+from . import utils
 
 try:
     from ._version import __version__
@@ -11,5 +12,6 @@
     "pqm_pvalue",
     "pqm_chi2",
     "test",
+    "utils",
     "__version__",
 )

src/pqm/utils.py

@@ -1,6 +1,6 @@
 import numpy as np
 import torch
-from scipy.stats import chi2
+from scipy.stats import chi2, beta
 from scipy.spatial.distance import cdist
 from tqdm import tqdm
 
@@ -13,6 +13,7 @@
     "_sample_reference_indices_torch",
     "_compute_counts_torch",
     "permute_test",
+    "max_test",
 )
 
 
@@ -272,3 +273,25 @@ def permute_test(f, x, y, n_permute=100, n_rerun=100, measure=np.mean):
         permute_stats.append(measure(list(f(x, y) for _ in range(n_rerun))))
 
     return test_stat, permute_stats
+
+
+def max_test(chi2_stats, dof):
+    """
+    Compute a p-value for the maximum of a set of chi2 values, under the
+    assumption they were drawn from a chi2_dof distribution.
+
+    Parameters
+    ----------
+    chi2_stats : np.ndarray
+        Test statistics.
+    dof : int
+        Degrees of freedom.
+
+    Returns
+    -------
+    float
+        p-value.
+    """
+    max_p_value = chi2.cdf(np.max(chi2_stats), dof)
+    p_value = 1 - (max_p_value) ** len(chi2_stats)
+    return p_value

tests/test_utils.py

@@ -0,0 +1,18 @@
+import numpy as np
+from scipy.stats import chi2
+
+from pqm import utils
+import pytest
+
+
+@pytest.mark.parametrize("dof", [2, 50, 100])
+def test_max_test(dof):
+    np.random.seed(42)
+    c = chi2.ppf(np.random.uniform(size=10), df=dof)
+
+    p = utils.max_test(c, dof)
+    assert p > 1e-4 and p < 0.9999, "p-value out of range, expected U(0,1)"
+
+    c[0] = 1000
+    p = utils.max_test(c, dof)
+    assert p < 1e-4, "p-value should be very small, expected ~0"