feat: Add toy calculator, empirical distribution, and toy example not…

…ebook (#790)

* Add EmpiricalDistribution and ToyCalculator classes
* Move top level infer functions under 'infer.calculators'
* Add tests for EmpiricalDistribution and ToyCalculator
* Add EmpiricalDistribution and ToyCalculator to docs
* Add example notebook on how to use the ToyCalculator

Co-authored-by: Giordon Stark <[email protected]>
Co-authored-by: Matthew Feickert <[email protected]>

docs/api.rst

@@ -128,10 +128,13 @@ Inference
    mle.fit
    mle.fixed_poi_fit
    hypotest
+   intervals.upperlimit
    calculators.generate_asimov_data
    calculators.AsymptoticTestStatDistribution
+   calculators.EmpiricalDistribution
    calculators.AsymptoticCalculator
-   intervals.upperlimit
+   calculators.ToyCalculator
+   utils.create_calculator
 
 Exceptions
 ----------

src/pyhf/cli/infer.py

@@ -138,6 +138,9 @@ def fit(
 @click.option('-p', '--patch', multiple=True)
 @click.option('--testpoi', default=1.0)
 @click.option('--teststat', type=click.Choice(['q', 'qtilde']), default='qtilde')
+@click.option(
+    '--calctype', type=click.Choice(['asymptotics', 'toybased']), default='asymptotics'
+)
 @click.option(
     '--backend',
     type=click.Choice(['numpy', 'pytorch', 'tensorflow', 'jax', 'np', 'torch', 'tf']),
@@ -160,6 +163,7 @@ def cls(
     teststat,
     backend,
     optimizer,
+    calctype,
     optconf,
 ):
     """
@@ -219,7 +223,12 @@ def cls(
         set_backend(tensorlib, new_optimizer(**optconf))
 
     result = hypotest(
-        testpoi, ws.data(model), model, qtilde=is_qtilde, return_expected_set=True
+        testpoi,
+        ws.data(model),
+        model,
+        qtilde=is_qtilde,
+        calctype=calctype,
+        return_expected_set=True,
     )
     result = {
         'CLs_obs': tensorlib.tolist(result[0]),

src/pyhf/infer/__init__.py

@@ -1,7 +1,7 @@
 """Inference for Statistical Models."""
 
+from . import utils
 from .. import get_backend
-from .calculators import AsymptoticCalculator
 
 
 def hypotest(
@@ -12,11 +12,17 @@ def hypotest(
     par_bounds=None,
     fixed_params=None,
     qtilde=True,
+    calctype="asymptotics",
+    return_tail_probs=False,
+    return_expected=False,
+    return_expected_set=False,
     **kwargs,
 ):
     r"""
     Compute :math:`p`-values and test statistics for a single value of the parameter of interest.
 
+    See :py:class:`~pyhf.infer.calculators.AsymptoticCalculator` and :py:class:`~pyhf.infer.calculators.ToyCalculator` on additional keyword arguments to be specified.
+
     Example:
         >>> import pyhf
         >>> pyhf.set_backend("numpy")
@@ -25,9 +31,9 @@ def hypotest(
         ... )
         >>> observations = [51, 48]
         >>> data = pyhf.tensorlib.astensor(observations + model.config.auxdata)
-        >>> test_poi = 1.0
+        >>> mu_test = 1.0
         >>> CLs_obs, CLs_exp_band = pyhf.infer.hypotest(
-        ...     test_poi, data, model, qtilde=True, return_expected_set=True
+        ...     mu_test, data, model, qtilde=True, return_expected_set=True
         ... )
         >>> CLs_obs
         array(0.05251497)
@@ -41,14 +47,13 @@ def hypotest(
         init_pars (:obj:`tensor`): The initial parameter values to be used for minimization
         par_bounds (:obj:`tensor`): The parameter value bounds to be used for minimization
         fixed_params (:obj:`tensor`): Whether to fix the parameter to the init_pars value during minimization
-        qtilde (:obj:`bool`): When ``True`` use :func:`~pyhf.infer.test_statistics.qmu_tilde`
-         as the test statistic.
-         When ``False`` use :func:`~pyhf.infer.test_statistics.qmu`.
-
-    Keyword Args:
-        return_tail_probs (bool): Bool for returning :math:`\mathrm{CL}_{s+b}` and :math:`\mathrm{CL}_{b}`
-        return_expected (bool): Bool for returning :math:`\mathrm{CL}_{\mathrm{exp}}`
-        return_expected_set (bool): Bool for returning the :math:`(-2,-1,0,1,2)\sigma` :math:`\mathrm{CL}_{\mathrm{exp}}` --- the "Brazil band"
+        qtilde (:obj:`bool`): When ``True`` perform the calculation using the alternative
+         test statistic, :math:`\tilde{q}_{\mu}`, as defined under the Wald
+         approximation in Equation (62) of :xref:`arXiv:1007.1727`.
+        calctype (:obj:`str`): The calculator to create. Choose either 'asymptotics' (default) or 'toybased'.
+        return_tail_probs (:obj:`bool`): Bool for returning :math:`\mathrm{CL}_{s+b}` and :math:`\mathrm{CL}_{b}`
+        return_expected (:obj:`bool`): Bool for returning :math:`\mathrm{CL}_{\mathrm{exp}}`
+        return_expected_set (:obj:`bool`): Bool for returning the :math:`(-2,-1,0,1,2)\sigma` :math:`\mathrm{CL}_{\mathrm{exp}}` --- the "Brazil band"
 
     Returns:
         Tuple of Floats and lists of Floats:
@@ -127,9 +132,17 @@ def hypotest(
     par_bounds = par_bounds or pdf.config.suggested_bounds()
     fixed_params = fixed_params or pdf.config.suggested_fixed()
 
-    calc = AsymptoticCalculator(
-        data, pdf, init_pars, par_bounds, fixed_params, qtilde=qtilde
+    calc = utils.create_calculator(
+        calctype,
+        data,
+        pdf,
+        init_pars,
+        par_bounds,
+        fixed_params,
+        qtilde=qtilde,
+        **kwargs,
     )
+
     teststat = calc.teststatistic(poi_test)
     sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test)
 
@@ -146,9 +159,9 @@ def hypotest(
     )
 
     _returns = [CLs]
-    if kwargs.get('return_tail_probs'):
+    if return_tail_probs:
         _returns.append([CLsb, CLb])
-    if kwargs.get('return_expected_set'):
+    if return_expected_set:
         CLs_exp = []
         for n_sigma in [2, 1, 0, -1, -2]:
 
@@ -158,10 +171,10 @@ def hypotest(
                 expected_bonly_teststat
             ) / b_only_distribution.pvalue(expected_bonly_teststat)
             CLs_exp.append(tensorlib.astensor(CLs))
-        if kwargs.get('return_expected'):
+        if return_expected:
             _returns.append(CLs_exp[2])
         _returns.append(CLs_exp)
-    elif kwargs.get('return_expected'):
+    elif return_expected:
         n_sigma = 0
         expected_bonly_teststat = b_only_distribution.expected_value(n_sigma)