feat: Add discovery p0 test statistic (#1232)

* Add discovery test statistic q0
* Teach the calculators to switch the asimov mu depending on the test statistic (to be refactored in later PR)
* Add tests for discovery test stat
* Add validation ROOT files and tests

Co-authored-by: Nikolai Hartmann <[email protected]>

docs/api.rst

@@ -122,6 +122,7 @@ Inference
 .. autosummary::
    :toctree: _generated/
 
+   test_statistics.q0
    test_statistics.qmu
    test_statistics.qmu_tilde
    test_statistics.tmu

src/pyhf/infer/__init__.py

@@ -153,18 +153,27 @@ def hypotest(
         tensorlib.astensor(CLs),
     )
 
-    _returns = [CLs]
+    is_q0 = kwargs.get('test_stat', 'qtilde') == 'q0'
+
+    _returns = [CLsb if is_q0 else CLs]
     if return_tail_probs:
-        _returns.append([CLsb, CLb])
+        if is_q0:
+            _returns.append([CLb])
+        else:
+            _returns.append([CLsb, CLb])
     if return_expected_set:
         CLs_exp = []
         for n_sigma in [2, 1, 0, -1, -2]:
 
             expected_bonly_teststat = b_only_distribution.expected_value(n_sigma)
 
-            CLs = sig_plus_bkg_distribution.pvalue(
-                expected_bonly_teststat
-            ) / b_only_distribution.pvalue(expected_bonly_teststat)
+            if is_q0:
+                # despite the name in this case this is the discovery p value
+                CLs = sig_plus_bkg_distribution.pvalue(expected_bonly_teststat)
+            else:
+                CLs = sig_plus_bkg_distribution.pvalue(
+                    expected_bonly_teststat
+                ) / b_only_distribution.pvalue(expected_bonly_teststat)
             CLs_exp.append(tensorlib.astensor(CLs))
         if return_expected:
             _returns.append(CLs_exp[2])
@@ -173,10 +182,16 @@ def hypotest(
         n_sigma = 0
         expected_bonly_teststat = b_only_distribution.expected_value(n_sigma)
 
-        CLs = sig_plus_bkg_distribution.pvalue(
-            expected_bonly_teststat
-        ) / b_only_distribution.pvalue(expected_bonly_teststat)
+        if is_q0:
+            # despite the name in this case this is the discovery p value
+            CLs = sig_plus_bkg_distribution.pvalue(expected_bonly_teststat)
+        else:
+            CLs = sig_plus_bkg_distribution.pvalue(
+                expected_bonly_teststat
+            ) / b_only_distribution.pvalue(expected_bonly_teststat)
+
         _returns.append(tensorlib.astensor(CLs))
+
     # Enforce a consistent return type of the observed CLs
     return tuple(_returns) if len(_returns) > 1 else _returns[0]
 

src/pyhf/infer/calculators.py

@@ -267,7 +267,8 @@ def teststatistic(self, poi_test):
         )
         sqrtqmu_v = tensorlib.sqrt(qmu_v)
 
-        asimov_mu = 0.0
+        asimov_mu = 1.0 if self.test_stat == 'q0' else 0.0
+
         asimov_data = generate_asimov_data(
             asimov_mu,
             self.data,
@@ -286,7 +287,7 @@ def teststatistic(self, poi_test):
         )
         self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)
 
-        if self.test_stat == "q":
+        if self.test_stat in ["q", "q0"]:  # qmu or q0
             teststat = sqrtqmu_v - self.sqrtqmuA_v
         else:  # qtilde
 
@@ -528,7 +529,7 @@ def distributions(self, poi_test, track_progress=None):
         signal_sample = signal_pdf.sample(sample_shape)
 
         bkg_pars = self.pdf.config.suggested_init()
-        bkg_pars[self.pdf.config.poi_index] = 0.0
+        bkg_pars[self.pdf.config.poi_index] = 1.0 if self.test_stat == 'q0' else 0.0
         bkg_pdf = self.pdf.make_pdf(tensorlib.astensor(bkg_pars))
         bkg_sample = bkg_pdf.sample(sample_shape)
 

src/pyhf/infer/test_statistics.py

@@ -290,3 +290,66 @@ def tmu_tilde(mu, data, pdf, init_pars, par_bounds, fixed_params):
             + 'Use the tmu test statistic (pyhf.infer.test_statistics.tmu) instead.'
         )
     return _tmu_like(mu, data, pdf, init_pars, par_bounds, fixed_params)
+
+
+def q0(mu, data, pdf, init_pars, par_bounds, fixed_params):
+    r"""
+    The test statistic, :math:`q_{0}`, for discovery of a positive signal
+    as defined in Equation (12) in :xref:`arXiv:1007.1727`, for :math:`\mu=0`.
+
+    .. math::
+       :nowrap:
+
+       \begin{equation}
+         q_{0} = \left\{\begin{array}{ll}
+         -2\ln\lambda\left(0\right), &\hat{\mu} \ge 0,\\
+         0, & \hat{\mu} < 0,
+         \end{array}\right.
+       \end{equation}
+
+    Example:
+        >>> import pyhf
+        >>> pyhf.set_backend("numpy")
+        >>> model = pyhf.simplemodels.hepdata_like(
+        ...     signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]
+        ... )
+        >>> observations = [60, 65]
+        >>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
+        >>> test_mu = 0.0
+        >>> init_pars = model.config.suggested_init()
+        >>> par_bounds = model.config.suggested_bounds()
+        >>> fixed_params = model.config.suggested_fixed()
+        >>> pyhf.infer.test_statistics.q0(test_mu, data, model, init_pars, par_bounds, fixed_params)
+        array(2.98339447)
+
+    Args:
+        mu (Number or Tensor): The signal strength parameter (must be set to zero)
+        data (Tensor): The data to be considered
+        pdf (~pyhf.pdf.Model): The HistFactory statistical model used in the likelihood ratio calculation
+        init_pars (:obj:`list`): Values to initialize the model parameters at for the fit
+        par_bounds (:obj:`list` of :obj:`list`\s or :obj:`tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
+        fixed_params (:obj:`list`): Parameters held constant in the fit
+
+    Returns:
+        Float: The calculated test statistic, :math:`q_{0}`
+    """
+
+    if pdf.config.poi_index is None:
+        raise UnspecifiedPOI(
+            'No POI is defined. A POI is required for profile likelihood based test statistics.'
+        )
+    if mu != 0.0:
+        log.warning(
+            'q0 test statistic only used for fit configuration with POI set to zero. Setting mu=0.'
+        )
+        mu = 0.0
+
+    tensorlib, optimizer = get_backend()
+
+    tmu_like_stat, (_, muhatbhat) = _tmu_like(
+        mu, data, pdf, init_pars, par_bounds, fixed_params, return_fitted_pars=True
+    )
+    q0_stat = tensorlib.where(
+        muhatbhat[pdf.config.poi_index] < 0, tensorlib.astensor(0.0), tmu_like_stat
+    )
+    return q0_stat

src/pyhf/infer/utils.py

@@ -2,7 +2,7 @@
 
 from .calculators import AsymptoticCalculator, ToyCalculator
 from ..exceptions import InvalidTestStatistic
-from .test_statistics import qmu, qmu_tilde
+from .test_statistics import q0, qmu, qmu_tilde
 
 import logging
 
@@ -49,12 +49,19 @@ def create_calculator(calctype, *args, **kwargs):
 
 def get_test_stat(name):
     """
-    Get the test statistic function, :func:`~pyhf.infer.test_statistics.qmu` or
-    :func:`~pyhf.infer.test_statistics.qmu_tilde`, by name.
+    Get the test statistic function by name. The following test statistics are supported:
+
+    - :func:`~pyhf.infer.test_statistics.q0`
+    - :func:`~pyhf.infer.test_statistics.qmu`
+    - :func:`~pyhf.infer.test_statistics.qmu_tilde`
 
     Example:
 
         >>> from pyhf.infer import utils, test_statistics
+        >>> utils.get_test_stat("q0")
+        <function q0 at 0x...>
+        >>> utils.get_test_stat("q0") == test_statistics.q0
+        True
         >>> utils.get_test_stat("q")
         <function qmu at 0x...>
         >>> utils.get_test_stat("q") == test_statistics.qmu
@@ -72,6 +79,7 @@ def get_test_stat(name):
         callable: The test statistic function
     """
     _mapping = {
+        "q0": q0,
         "q": qmu,
         "qtilde": qmu_tilde,
     }

tests/test_infer.py

@@ -93,40 +93,47 @@ def test_hypotest_poi_outofbounds(tmpdir, hypotest_args):
         pyhf.infer.hypotest(10.1, data, pdf)
 
 
-def test_hypotest_return_tail_probs(tmpdir, hypotest_args):
+@pytest.mark.parametrize('test_stat', ['q0', 'q', 'qtilde'])
+def test_hypotest_return_tail_probs(tmpdir, hypotest_args, test_stat):
     """
     Check that the return structure of pyhf.infer.hypotest with the
     return_tail_probs keyword arg is as expected
     """
     tb = pyhf.tensorlib
 
-    kwargs = {'return_tail_probs': True}
+    kwargs = {'return_tail_probs': True, 'test_stat': test_stat}
     result = pyhf.infer.hypotest(*hypotest_args, **kwargs)
     # CLs_obs, [CL_sb, CL_b]
     assert len(list(result)) == 2
     assert isinstance(result[0], type(tb.astensor(result[0])))
-    assert len(result[1]) == 2
+    assert len(result[1]) == 1 if test_stat == 'q0' else 2
     assert check_uniform_type(result[1])
 
 
-def test_hypotest_return_expected(tmpdir, hypotest_args):
+@pytest.mark.parametrize('test_stat', ['q0', 'q', 'qtilde'])
+def test_hypotest_return_expected(tmpdir, hypotest_args, test_stat):
     """
     Check that the return structure of pyhf.infer.hypotest with the
     additon of the return_expected keyword arg is as expected
     """
     tb = pyhf.tensorlib
 
-    kwargs = {'return_tail_probs': True, 'return_expected': True}
+    kwargs = {
+        'return_tail_probs': True,
+        'return_expected': True,
+        'test_stat': test_stat,
+    }
     result = pyhf.infer.hypotest(*hypotest_args, **kwargs)
     # CLs_obs, [CLsb, CLb], CLs_exp
     assert len(list(result)) == 3
     assert isinstance(result[0], type(tb.astensor(result[0])))
-    assert len(result[1]) == 2
+    assert len(result[1]) == 1 if test_stat == 'q0' else 2
     assert check_uniform_type(result[1])
     assert isinstance(result[2], type(tb.astensor(result[2])))
 
 
-def test_hypotest_return_expected_set(tmpdir, hypotest_args):
+@pytest.mark.parametrize('test_stat', ['q0', 'q', 'qtilde'])
+def test_hypotest_return_expected_set(tmpdir, hypotest_args, test_stat):
     """
     Check that the return structure of pyhf.infer.hypotest with the
     additon of the return_expected_set keyword arg is as expected
@@ -137,12 +144,13 @@ def test_hypotest_return_expected_set(tmpdir, hypotest_args):
         'return_tail_probs': True,
         'return_expected': True,
         'return_expected_set': True,
+        'test_stat': test_stat,
     }
     result = pyhf.infer.hypotest(*hypotest_args, **kwargs)
     # CLs_obs, [CLsb, CLb], CLs_exp, CLs_exp @[-2, -1, 0, +1, +2]sigma
     assert len(list(result)) == 4
     assert isinstance(result[0], type(tb.astensor(result[0])))
-    assert len(result[1]) == 2
+    assert len(result[1]) == 1 if test_stat == 'q0' else 2
     assert check_uniform_type(result[1])
     assert isinstance(result[2], type(tb.astensor(result[2])))
     assert len(result[3]) == 5

tests/test_teststats.py

@@ -4,6 +4,25 @@
 import logging
 
 
+def test_q0(caplog):
+    mu = 1.0
+    model = pyhf.simplemodels.hepdata_like([6], [9], [3])
+    data = [9] + model.config.auxdata
+    init_pars = model.config.suggested_init()
+    par_bounds = model.config.suggested_bounds()
+    fixed_params = model.config.suggested_fixed()
+
+    with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"):
+        pyhf.infer.test_statistics.q0(
+            mu, data, model, init_pars, par_bounds, fixed_params
+        )
+        assert (
+            "q0 test statistic only used for fit configuration with POI set to zero"
+            in caplog.text
+        )
+        caplog.clear()
+
+
 def test_qmu(caplog):
     mu = 1.0
     model = pyhf.simplemodels.hepdata_like([6], [9], [3])
@@ -99,6 +118,15 @@ def test_no_poi_test_stats():
     par_bounds = model.config.suggested_bounds()
     fixed_params = model.config.suggested_fixed()
 
+    with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo:
+        pyhf.infer.test_statistics.q0(
+            test_poi, data, model, init_pars, par_bounds, fixed_params
+        )
+    assert (
+        "No POI is defined. A POI is required for profile likelihood based test statistics."
+        in str(excinfo.value)
+    )
+
     with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo:
         pyhf.infer.test_statistics.qmu(
             test_poi, data, model, init_pars, par_bounds, fixed_params