Merge pull request #10 from Ciela-Institute/fixnumericalstability

fix numerical stability for chi2 rescaling

.gitignore

@@ -1,2 +1,5 @@
 **__pycache__
-**_version.py
+**_version.py
+**.gz
+**.ipynb_checkpoints
+**ubyte

src/pqm/pqm.py

@@ -135,7 +135,10 @@ def pqm_chi2(
     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)
+        if chi2_stat / dof < 10:
+            cp = chi2.sf(chi2_stat, dof)
+            chi2_stat = chi2.isf(cp, num_refs - 1)
+        else:
+            chi2_stat = chi2_stat * (num_refs - 1) / dof
         dof = num_refs - 1
     return chi2_stat, dof

tests/test_gaussian.py

@@ -43,4 +43,4 @@ def test_fail_chi2():
 
         new.append(pqm_chi2(x_samples, y_samples, num_refs=100))
     new = np.array(new)
-    assert np.mean(new[:, 0]) / 99 > 10
+    assert np.mean(new[:, 0]) / 99 > 2