diff --git a/src/pqm/pqm.py b/src/pqm/pqm.py index 3a7cfa1..a5bbdfa 100644 --- a/src/pqm/pqm.py +++ b/src/pqm/pqm.py @@ -1,7 +1,7 @@ from typing import Optional import numpy as np -from scipy.stats import chi2_contingency +from scipy.stats import chi2_contingency, chi2 from scipy.spatial import KDTree __all__ = ("pqm_chi2", "pqm_pvalue") @@ -133,4 +133,9 @@ def pqm_chi2( for _ in range(bootstrap) ] chi2_stat, _, dof, _ = _pqm_test(x_samples, y_samples, num_refs, whiten) + if dof != num_refs - 1: + # Rescale chi2 to new value which has the same cumulative probability + cp = chi2.cdf(chi2_stat, dof) + chi2_stat = chi2.ppf(cp, num_refs - 1) + dof = num_refs - 1 return chi2_stat, dof diff --git a/tests/test_gaussian.py b/tests/test_gaussian.py index 2c113b8..4e3eb1c 100644 --- a/tests/test_gaussian.py +++ b/tests/test_gaussian.py @@ -1,8 +1,8 @@ import numpy as np -from pqm import pqm_pvalue +from pqm import pqm_pvalue, pqm_chi2 -def test(): +def test_pass_pvalue(): new = [] for _ in range(100): y_samples = np.random.normal(size=(500, 50)) @@ -11,3 +11,36 @@ def test(): new.append(pqm_pvalue(x_samples, y_samples)) assert np.abs(np.mean(new) - 0.5) < 0.15 + + +def test_pass_chi2(): + new = [] + for _ in range(100): + y_samples = np.random.normal(size=(500, 50)) + x_samples = np.random.normal(size=(250, 50)) + + new.append(pqm_chi2(x_samples, y_samples, num_refs=100)) + new = np.array(new) + assert np.abs(np.mean(new[:, 0]) / 99 - 1) < 0.15 + + +def test_fail_pvalue(): + new = [] + for _ in range(100): + y_samples = np.random.normal(size=(500, 50)) + x_samples = np.random.normal(size=(250, 50)) + 0.5 + + new.append(pqm_pvalue(x_samples, y_samples)) + + assert np.mean(new) < 1e-3 + + +def test_fail_chi2(): + new = [] + for _ in range(100): + y_samples = np.random.normal(size=(500, 50)) + x_samples = np.random.normal(size=(250, 50)) + 0.5 + + new.append(pqm_chi2(x_samples, y_samples, num_refs=100)) + new = np.array(new) + assert np.mean(new[:, 0]) / 99 > 10