refactor: Factor out tmu-like and qmu-like test statistics (#1005)

* Add _tmu_like and _qmu_like functions to separate behavior
* Add qmu, qmu_tilde, tmu, and tmu_tilde functions for more explicit API
* Add functions to the docs
* Add tests for new API and warnings raised

Author: lukasheinrich

docs/api.rst

@@ -122,6 +122,9 @@ Inference
 
    hypotest
    test_statistics.qmu
+   test_statistics.qmu_tilde
+   test_statistics.tmu
+   test_statistics.tmu_tilde
    mle.twice_nll
    mle.fit
    mle.fixed_poi_fit

src/pyhf/infer/calculators.py

@@ -9,7 +9,7 @@
 """
 from .mle import fixed_poi_fit
 from .. import get_backend
-from .test_statistics import qmu
+from .test_statistics import qmu, qmu_tilde
 
 
 def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
@@ -155,14 +155,21 @@ def teststatistic(self, poi_test):
 
         """
         tensorlib, _ = get_backend()
-        qmu_v = qmu(poi_test, self.data, self.pdf, self.init_pars, self.par_bounds)
+
+        teststat_func = qmu_tilde if self.qtilde else qmu
+
+        qmu_v = teststat_func(
+            poi_test, self.data, self.pdf, self.init_pars, self.par_bounds
+        )
         sqrtqmu_v = tensorlib.sqrt(qmu_v)
 
         asimov_mu = 0.0
         asimov_data = generate_asimov_data(
             asimov_mu, self.data, self.pdf, self.init_pars, self.par_bounds
         )
-        qmuA_v = qmu(poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds)
+        qmuA_v = teststat_func(
+            poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds
+        )
         self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)
 
         if not self.qtilde:  # qmu

src/pyhf/infer/test_statistics.py

@@ -2,6 +2,51 @@
 from .mle import fixed_poi_fit, fit
 from ..exceptions import UnspecifiedPOI
 
+import logging
+
+log = logging.getLogger(__name__)
+
+
+def _qmu_like(mu, data, pdf, init_pars, par_bounds):
+    """
+    Clipped version of _tmu_like where the returned test statistic
+    is 0 if muhat > 0 else tmu_like_stat.
+
+    If the lower bound of the POI is 0 this automatically implments
+    qmu_tilde. Otherwise this is qmu (no tilde).
+    """
+    tensorlib, optimizer = get_backend()
+    tmu_like_stat, (_, muhatbhat) = _tmu_like(
+        mu, data, pdf, init_pars, par_bounds, return_fitted_pars=True
+    )
+    qmu_like_stat = tensorlib.where(
+        muhatbhat[pdf.config.poi_index] > mu, tensorlib.astensor(0.0), tmu_like_stat
+    )
+    return qmu_like_stat
+
+
+def _tmu_like(mu, data, pdf, init_pars, par_bounds, return_fitted_pars=False):
+    """
+    Basic Profile Likelihood test statistic.
+
+    If the lower bound of the POI is 0 this automatically implments
+    tmu_tilde. Otherwise this is tmu (no tilde).
+    """
+    tensorlib, optimizer = get_backend()
+    mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit(
+        mu, data, pdf, init_pars, par_bounds, return_fitted_val=True
+    )
+    muhatbhat, unconstrained_fit_lhood_val = fit(
+        data, pdf, init_pars, par_bounds, return_fitted_val=True
+    )
+    log_likelihood_ratio = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val
+    tmu_like_stat = tensorlib.astensor(
+        tensorlib.clip(log_likelihood_ratio, 0.0, max_value=None)
+    )
+    if return_fitted_pars:
+        return tmu_like_stat, (mubhathat, muhatbhat)
+    return tmu_like_stat
+
 
 def qmu(mu, data, pdf, init_pars, par_bounds):
     r"""
@@ -30,8 +75,9 @@ def qmu(mu, data, pdf, init_pars, par_bounds):
         >>> test_mu = 1.0
         >>> init_pars = model.config.suggested_init()
         >>> par_bounds = model.config.suggested_bounds()
+        >>> par_bounds[model.config.poi_index] = [-10.0, 10.0]
         >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds)
-        3.938244920380498
+        array(3.9549891)
 
     Args:
         mu (Number or Tensor): The signal strength parameter
@@ -47,16 +93,136 @@ def qmu(mu, data, pdf, init_pars, par_bounds):
         raise UnspecifiedPOI(
             'No POI is defined. A POI is required for profile likelihood based test statistics.'
         )
+    if par_bounds[pdf.config.poi_index][0] == 0:
+        log.warning(
+            'qmu test statistic used for fit configuration with POI bounded at zero.\n'
+            + 'Use the qmu_tilde test statistic (pyhf.infer.test_statistics.qmu_tilde) instead.'
+        )
+    return _qmu_like(mu, data, pdf, init_pars, par_bounds)
 
-    tensorlib, optimizer = get_backend()
-    mubhathat, fixed_poi_fit_lhood_val = fixed_poi_fit(
-        mu, data, pdf, init_pars, par_bounds, return_fitted_val=True
-    )
-    muhatbhat, unconstrained_fit_lhood_val = fit(
-        data, pdf, init_pars, par_bounds, return_fitted_val=True
-    )
-    qmu = fixed_poi_fit_lhood_val - unconstrained_fit_lhood_val
-    qmu = tensorlib.where(
-        muhatbhat[pdf.config.poi_index] > mu, tensorlib.astensor(0.0), qmu
-    )
-    return tensorlib.clip(qmu, 0, max_value=None)
+
+def qmu_tilde(mu, data, pdf, init_pars, par_bounds):
+    r"""
+    The test statistic, :math:`\tilde{q}_{\mu}`, for establishing an upper
+    limit on the strength parameter, :math:`\mu`, for models with
+    bounded POI, as defiend in Equation (16) in :xref:`arXiv:1007.1727`.
+
+    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 = [51, 48]
+        >>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
+        >>> test_mu = 1.0
+        >>> init_pars = model.config.suggested_init()
+        >>> par_bounds = model.config.suggested_bounds()
+        >>> pyhf.infer.test_statistics.qmu_tilde(test_mu, data, model, init_pars, par_bounds)
+        array(3.93824492)
+
+    Args:
+        mu (Number or Tensor): The signal strength parameter
+        data (Tensor): The data to be considered
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+        init_pars (`list`): Values to initialize the model parameters at for the fit
+        par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
+
+    Returns:
+        Float: The calculated test statistic, :math:`\tilde{q}_{\mu}`
+    """
+    if pdf.config.poi_index is None:
+        raise UnspecifiedPOI(
+            'No POI is defined. A POI is required for profile likelihood based test statistics.'
+        )
+    if par_bounds[pdf.config.poi_index][0] != 0:
+        log.warning(
+            'qmu_tilde test statistic used for fit configuration with POI not bounded at zero.\n'
+            + 'Use the qmu test statistic (pyhf.infer.test_statistics.qmu) instead.'
+        )
+    return _qmu_like(mu, data, pdf, init_pars, par_bounds)
+
+
+def tmu(mu, data, pdf, init_pars, par_bounds):
+    r"""
+    The test statistic, :math:`t_{\mu}`, for establishing a two-sided
+    interval on the strength parameter, :math:`\mu`, as defiend in Equation (10)
+    in :xref:`arXiv:1007.1727`.
+
+    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 = [51, 48]
+        >>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
+        >>> test_mu = 1.0
+        >>> init_pars = model.config.suggested_init()
+        >>> par_bounds = model.config.suggested_bounds()
+        >>> par_bounds[model.config.poi_index] = [-10.0, 10.0]
+        >>> pyhf.infer.test_statistics.tmu(test_mu, data, model, init_pars, par_bounds)
+        array(3.9549891)
+
+    Args:
+        mu (Number or Tensor): The signal strength parameter
+        data (Tensor): The data to be considered
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+        init_pars (`list`): Values to initialize the model parameters at for the fit
+        par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
+
+    Returns:
+        Float: The calculated test statistic, :math:`t_{\mu}`
+    """
+    if pdf.config.poi_index is None:
+        raise UnspecifiedPOI(
+            'No POI is defined. A POI is required for profile likelihood based test statistics.'
+        )
+    if par_bounds[pdf.config.poi_index][0] == 0:
+        log.warning(
+            'tmu test statistic used for fit configuration with POI bounded at zero.\n'
+            + 'Use the tmu_tilde test statistic (pyhf.infer.test_statistics.tmu_tilde) instead.'
+        )
+    return _tmu_like(mu, data, pdf, init_pars, par_bounds)
+
+
+def tmu_tilde(mu, data, pdf, init_pars, par_bounds):
+    r"""
+    The test statistic, :math:`t_{\mu}`, for establishing a two-sided
+    interval on the strength parameter, :math:`\mu`, for models with
+    bounded POI, as defiend in Equation (11) in :xref:`arXiv:1007.1727`.
+
+    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 = [51, 48]
+        >>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
+        >>> test_mu = 1.0
+        >>> init_pars = model.config.suggested_init()
+        >>> par_bounds = model.config.suggested_bounds()
+        >>> pyhf.infer.test_statistics.tmu_tilde(test_mu, data, model, init_pars, par_bounds)
+        array(3.93824492)
+
+    Args:
+        mu (Number or Tensor): The signal strength parameter
+        data (Tensor): The data to be considered
+        pdf (~pyhf.pdf.Model): The statistical model adhering to the schema model.json
+        init_pars (`list`): Values to initialize the model parameters at for the fit
+        par_bounds (`list` of `list`\s or `tuple`\s): The extrema of values the model parameters are allowed to reach in the fit
+
+    Returns:
+        Float: The calculated test statistic, :math:`\tilde{t}_{\mu}`
+    """
+    if pdf.config.poi_index is None:
+        raise UnspecifiedPOI(
+            'No POI is defined. A POI is required for profile likelihood based test statistics.'
+        )
+    if par_bounds[pdf.config.poi_index][0] != 0:
+        log.warning(
+            'tmu_tilde test statistic used for fit configuration with POI not bounded at zero.\n'
+            + 'Use the tmu test statistic (pyhf.infer.test_statistics.tmu) instead.'
+        )
+    return _tmu_like(mu, data, pdf, init_pars, par_bounds)

tests/test_backend_consistency.py

@@ -58,7 +58,7 @@ def generate_source_poisson(n_bins):
 
 @pytest.mark.parametrize('n_bins', bins, ids=bin_ids)
 @pytest.mark.parametrize('invert_order', [False, True], ids=['normal', 'inverted'])
-def test_hypotest_q_mu(
+def test_hypotest_qmu_tilde(
     n_bins, invert_order, tolerance={'numpy': 1e-02, 'tensors': 5e-03}
 ):
     """
@@ -114,10 +114,10 @@ def test_hypotest_q_mu(
     for backend in backends:
         pyhf.set_backend(backend)
 
-        q_mu = pyhf.infer.test_statistics.qmu(
+        qmu_tilde = pyhf.infer.test_statistics.qmu_tilde(
             1.0, data, pdf, pdf.config.suggested_init(), pdf.config.suggested_bounds(),
         )
-        test_statistic.append(q_mu)
+        test_statistic.append(qmu_tilde)
 
     # compare to NumPy/SciPy
     test_statistic = np.array(test_statistic)

tests/test_teststats.py

@@ -0,0 +1,118 @@
+import pytest
+import pyhf
+import pyhf.infer.test_statistics
+import logging
+
+
+def test_qmu(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()
+
+    with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"):
+        pyhf.infer.test_statistics.qmu(mu, data, model, init_pars, par_bounds)
+        assert "WARNING  qmu test statistic used for fit" in caplog.text
+        caplog.clear()
+
+
+def test_qmu_tilde(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()
+
+    par_bounds[model.config.poi_index] = [-10, 10]
+    with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"):
+        pyhf.infer.test_statistics.qmu_tilde(mu, data, model, init_pars, par_bounds)
+        assert "WARNING  qmu_tilde test statistic used for fit" in caplog.text
+        caplog.clear()
+
+
+def test_tmu(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()
+    with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"):
+        pyhf.infer.test_statistics.tmu(mu, data, model, init_pars, par_bounds)
+        assert "WARNING  tmu test statistic used for fit" in caplog.text
+        caplog.clear()
+
+
+def test_tmu_tilde(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()
+
+    par_bounds[model.config.poi_index] = [-10, 10]
+    with caplog.at_level(logging.WARNING, "pyhf.infer.test_statistics"):
+        pyhf.infer.test_statistics.tmu_tilde(mu, data, model, init_pars, par_bounds)
+        assert "WARNING  tmu_tilde test statistic used for fit" in caplog.text
+        caplog.clear()
+
+
+def test_no_poi_test_stats():
+    spec = {
+        "channels": [
+            {
+                "name": "channel",
+                "samples": [
+                    {
+                        "name": "sample",
+                        "data": [10.0],
+                        "modifiers": [
+                            {
+                                "type": "normsys",
+                                "name": "shape",
+                                "data": {"hi": 0.5, "lo": 1.5},
+                            }
+                        ],
+                    },
+                ],
+            }
+        ]
+    }
+    model = pyhf.Model(spec, poi_name=None)
+
+    test_poi = 1.0
+    data = [12] + model.config.auxdata
+    init_pars = model.config.suggested_init()
+    par_bounds = model.config.suggested_bounds()
+
+    with pytest.raises(pyhf.exceptions.UnspecifiedPOI) as excinfo:
+        pyhf.infer.test_statistics.qmu(test_poi, data, model, init_pars, par_bounds)
+    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_tilde(
+            test_poi, data, model, init_pars, par_bounds
+        )
+    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.tmu(test_poi, data, model, init_pars, par_bounds)
+    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.tmu_tilde(
+            test_poi, data, model, init_pars, par_bounds
+        )
+    assert (
+        "No POI is defined. A POI is required for profile likelihood based test statistics."
+        in str(excinfo.value)
+    )