From ee183ebfa9c01bfeb3b3f06d82e5b4e9bf518753 Mon Sep 17 00:00:00 2001
From: Matthew Feickert <matthew.feickert@cern.ch>
Date: Fri, 9 Sep 2022 17:46:48 +0200
Subject: [PATCH] fix: Guard against nan in test stat calculation (#1993)

* Guard against nan from division by zero in pyhf.infer.calculators.AsymptoticCalculator.teststatistic.
* Add use of 'less than or equal to' to docs for tests stats to match equations 14
  and 16 of https://arxiv.org/abs/1007.1727.
* Add tests to ensure that nan conditions in Issue #529 and Issue #1992 are not possible.
---
 src/pyhf/infer/calculators.py     |  3 ++-
 src/pyhf/infer/test_statistics.py |  4 ++--
 tests/test_infer.py               | 34 +++++++++++++++++++++++++++++++
 3 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py
index 8e06a613cb..ed9642e085 100644
--- a/src/pyhf/infer/calculators.py
+++ b/src/pyhf/infer/calculators.py
@@ -418,8 +418,9 @@ def _false_case():
                 teststat = (qmu - qmu_A) / (2 * self.sqrtqmuA_v)
                 return teststat
 
+            # Use '<=' rather than '<' to avoid Issue #1992
             teststat = tensorlib.conditional(
-                (sqrtqmu_v < self.sqrtqmuA_v), _true_case, _false_case
+                (sqrtqmu_v <= self.sqrtqmuA_v), _true_case, _false_case
             )
         return tensorlib.astensor(teststat)
 
diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py
index e23e87fa75..3af773d42c 100644
--- a/src/pyhf/infer/test_statistics.py
+++ b/src/pyhf/infer/test_statistics.py
@@ -71,7 +71,7 @@ def qmu(mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=F
 
        \begin{equation}
           q_{\mu} = \left\{\begin{array}{ll}
-          -2\ln\lambda\left(\mu\right), &\hat{\mu} < \mu,\\
+          -2\ln\lambda\left(\mu\right), &\hat{\mu} \leq \mu,\\
           0, & \hat{\mu} > \mu
           \end{array}\right.
         \end{equation}
@@ -160,7 +160,7 @@ def qmu_tilde(
 
        \begin{equation}
           \tilde{q}_{\mu} = \left\{\begin{array}{ll}
-          -2\ln\tilde{\lambda}\left(\mu\right), &\hat{\mu} < \mu,\\
+          -2\ln\tilde{\lambda}\left(\mu\right), &\hat{\mu} \leq \mu,\\
           0, & \hat{\mu} > \mu
           \end{array}\right.
         \end{equation}
diff --git a/tests/test_infer.py b/tests/test_infer.py
index 92a3467585..8d06f57694 100644
--- a/tests/test_infer.py
+++ b/tests/test_infer.py
@@ -476,3 +476,37 @@ def test_fixed_poi(tmpdir, hypotest_args):
     pdf.config.param_set('mu').suggested_fixed = [True]
     with pytest.raises(pyhf.exceptions.InvalidModel):
         pyhf.infer.hypotest(*hypotest_args)
+
+
+def test_teststat_nan_guard():
+    # Example from Issue #1992
+    model = pyhf.simplemodels.uncorrelated_background(
+        signal=[1.0], bkg=[1.0], bkg_uncertainty=[1.0]
+    )
+    observations = [2]
+    test_poi = 0.0
+    data = observations + model.config.auxdata
+    init_pars = model.config.suggested_init()
+    par_bounds = model.config.suggested_bounds()
+    fixed_params = model.config.suggested_fixed()
+
+    test_stat = pyhf.infer.test_statistics.qmu_tilde(
+        test_poi, data, model, init_pars, par_bounds, fixed_params
+    )
+    assert test_stat == pytest.approx(0.0)
+    asymptotic_calculator = pyhf.infer.calculators.AsymptoticCalculator(
+        data, model, test_stat="qtilde"
+    )
+    # ensure not nan
+    assert ~np.isnan(asymptotic_calculator.teststatistic(test_poi))
+    assert asymptotic_calculator.teststatistic(test_poi) == pytest.approx(0.0)
+
+    # Example from Issue #529
+    model = pyhf.simplemodels.uncorrelated_background([0.005], [28.0], [5.0])
+    test_poi = 1.0
+    data = [28.0] + model.config.auxdata
+
+    test_results = pyhf.infer.hypotest(
+        test_poi, data, model, test_stat="qtilde", return_expected=True
+    )
+    assert all(~np.isnan(result) for result in test_results)