feat: Support configuring fixed values automatically from Model (#1051)

* Add support for automatically passing through fixed parameters from the model configuration into all of the minimization/inference API * Add tests for fixed values Co-authored-by: Carlo Gottardo <[email protected]>
scikit-hep · Sep 7, 2020 · 2cabb09 · 2cabb09
1 parent 516ee1b
commit 2cabb09
Show file tree

Hide file tree

Showing 10 changed files with 163 additions and 57 deletions.
diff --git a/docs/examples/notebooks/ImpactPlot.ipynb b/docs/examples/notebooks/ImpactPlot.ipynb
@@ -82,8 +82,17 @@
     "    _, model, data = make_model([\"CRtt_meff\"])\n",
     "\n",
     "    pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n",
+    "    \n",
+    "    constraints = constraints or []\n",
+    "    init_pars = model.config.suggested_init()\n",
+    "    fixed_params = model.config.suggested_fixed()\n",
+    "    \n",
+    "    for idx,fixed_val in constraints:\n",
+    "        init_pars[idx] = fixed_val\n",
+    "        fixed_params[idx] = True\n",
+    "    \n",
     "    result = pyhf.infer.mle.fit(\n",
-    "        data, model, fixed_vals=constraints, return_uncertainties=True\n",
+    "        data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n",
     "    )\n",
     "    bestfit = result[:, 0]\n",
     "    errors = result[:, 1]\n",

diff --git a/docs/examples/notebooks/pullplot.ipynb b/docs/examples/notebooks/pullplot.ipynb
@@ -72,8 +72,17 @@
     "    _, model, data = make_model([\"CRtt_meff\"])\n",
     "\n",
     "    pyhf.set_backend(\"numpy\", pyhf.optimize.minuit_optimizer(verbose=True))\n",
+    "    \n",
+    "    constraints = constraints or []\n",
+    "    init_pars = model.config.suggested_init()\n",
+    "    fixed_params = model.config.suggested_fixed()\n",
+    "    \n",
+    "    for idx,fixed_val in constraints:\n",
+    "        init_pars[idx] = fixed_val\n",
+    "        fixed_params[idx] = True\n",
+    "    \n",
     "    result = pyhf.infer.mle.fit(\n",
-    "        data, model, fixed_vals=constraints, return_uncertainties=True\n",
+    "        data, model, init_pars=init_pars, fixed_params=fixed_params, return_uncertainties=True\n",
     "    )\n",
     "    bestfit = result[:, 0]\n",
     "    errors = result[:, 1]\n",

diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py
@@ -6,7 +6,14 @@
 
 
 def hypotest(
-    poi_test, data, pdf, init_pars=None, par_bounds=None, qtilde=False, **kwargs
+    poi_test,
+    data,
+    pdf,
+    init_pars=None,
+    par_bounds=None,
+    fixed_params=None,
+    qtilde=False,
+    **kwargs,
 ):
     r"""
     Compute :math:`p`-values and test statistics for a single value of the parameter of interest.
@@ -32,8 +39,9 @@ def hypotest(
         poi_test (Number or Tensor): The value of the parameter of interest (POI)
         data (Number or Tensor): The data considered
         pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``
-        init_pars (Array or Tensor): The initial parameter values to be used for minimization
-        par_bounds (Array or Tensor): The parameter value bounds to be used for minimization
+        init_pars (`tensor`): The initial parameter values to be used for minimization
+        par_bounds (`tensor`): The parameter value bounds to be used for minimization
+        fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization
         qtilde (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`.
@@ -118,15 +126,19 @@ def hypotest(
     """
     init_pars = init_pars or pdf.config.suggested_init()
     par_bounds = par_bounds or pdf.config.suggested_bounds()
-    tensorlib, _ = get_backend()
+    fixed_params = fixed_params or pdf.config.suggested_fixed()
 
-    calc = AsymptoticCalculator(data, pdf, init_pars, par_bounds, qtilde=qtilde)
+    calc = AsymptoticCalculator(
+        data, pdf, init_pars, par_bounds, fixed_params, qtilde=qtilde
+    )
     teststat = calc.teststatistic(poi_test)
     sig_plus_bkg_distribution, b_only_distribution = calc.distributions(poi_test)
 
     CLsb = sig_plus_bkg_distribution.pvalue(teststat)
     CLb = b_only_distribution.pvalue(teststat)
     CLs = CLsb / CLb
+
+    tensorlib, _ = get_backend()
     # Ensure that all CL values are 0-d tensors
     CLsb, CLb, CLs = (
         tensorlib.astensor(CLsb),

diff --git a/src/pyhf/infer/calculators.py b/src/pyhf/infer/calculators.py
@@ -12,7 +12,7 @@
 from .test_statistics import qmu, qmu_tilde
 
 
-def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
+def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds, fixed_params):
     """
     Compute Asimov Dataset (expected yields at best-fit values) for a given POI value.
 
@@ -22,12 +22,15 @@ def generate_asimov_data(asimov_mu, data, pdf, init_pars, par_bounds):
         pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
         init_pars (`tensor`): The initial parameter values to be used for fitting.
         par_bounds (`tensor`): The parameter value bounds to be used for fitting.
+        fixed_params (`tensor`): Parameters to be held constant in the fit.
 
     Returns:
         Tensor: The Asimov dataset.
 
     """
-    bestfit_nuisance_asimov = fixed_poi_fit(asimov_mu, data, pdf, init_pars, par_bounds)
+    bestfit_nuisance_asimov = fixed_poi_fit(
+        asimov_mu, data, pdf, init_pars, par_bounds, fixed_params
+    )
     return pdf.expected_data(bestfit_nuisance_asimov)
 
 
@@ -118,7 +121,15 @@ def expected_value(self, nsigma):
 class AsymptoticCalculator(object):
     """The Asymptotic Calculator."""
 
-    def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
+    def __init__(
+        self,
+        data,
+        pdf,
+        init_pars=None,
+        par_bounds=None,
+        fixed_params=None,
+        qtilde=False,
+    ):
         """
         Asymptotic Calculator.
 
@@ -127,6 +138,7 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
             pdf (~pyhf.pdf.Model): The statistical model adhering to the schema ``model.json``.
             init_pars (`tensor`): The initial parameter values to be used for fitting.
             par_bounds (`tensor`): The parameter value bounds to be used for fitting.
+            fixed_params (`tensor`): Whether to fix the parameter to the init_pars value during minimization
             qtilde (`bool`): Whether to use qtilde as the test statistic.
 
         Returns:
@@ -137,6 +149,8 @@ def __init__(self, data, pdf, init_pars=None, par_bounds=None, qtilde=False):
         self.pdf = pdf
         self.init_pars = init_pars or pdf.config.suggested_init()
         self.par_bounds = par_bounds or pdf.config.suggested_bounds()
+        self.fixed_params = fixed_params or pdf.config.suggested_fixed()
+
         self.qtilde = qtilde
         self.sqrtqmuA_v = None
 
@@ -173,16 +187,31 @@ def teststatistic(self, poi_test):
         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
+            poi_test,
+            self.data,
+            self.pdf,
+            self.init_pars,
+            self.par_bounds,
+            self.fixed_params,
         )
         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
+            asimov_mu,
+            self.data,
+            self.pdf,
+            self.init_pars,
+            self.par_bounds,
+            self.fixed_params,
         )
         qmuA_v = teststat_func(
-            poi_test, asimov_data, self.pdf, self.init_pars, self.par_bounds
+            poi_test,
+            asimov_data,
+            self.pdf,
+            self.init_pars,
+            self.par_bounds,
+            self.fixed_params,
         )
         self.sqrtqmuA_v = tensorlib.sqrt(qmuA_v)
 

diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py
@@ -47,9 +47,9 @@ def twice_nll(pars, data, pdf):
     return -2 * pdf.logpdf(pars, data)
 
 
-def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
+def fit(data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs):
     r"""
-    Run a unconstrained maximum likelihood fit.
+    Run a maximum likelihood fit.
     This is done by minimizing the objective function :func:`~pyhf.infer.mle.twice_nll`
     of the model parameters given the observed data.
     This is used to produce the maximal likelihood :math:`L\left(\hat{\mu}, \hat{\boldsymbol{\theta}}\right)`
@@ -87,6 +87,7 @@ def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
         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
+        fixed_params (`list`): Parameters to be held constant in the fit.
         kwargs: Keyword arguments passed through to the optimizer API
 
     Returns:
@@ -96,10 +97,23 @@ def fit(data, pdf, init_pars=None, par_bounds=None, **kwargs):
     _, opt = get_backend()
     init_pars = init_pars or pdf.config.suggested_init()
     par_bounds = par_bounds or pdf.config.suggested_bounds()
-    return opt.minimize(twice_nll, data, pdf, init_pars, par_bounds, **kwargs)
+    fixed_params = fixed_params or pdf.config.suggested_fixed()
 
+    # get fixed vals from the model
+    fixed_vals = [
+        (index, init)
+        for index, (init, is_fixed) in enumerate(zip(init_pars, fixed_params))
+        if is_fixed
+    ]
 
-def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs):
+    return opt.minimize(
+        twice_nll, data, pdf, init_pars, par_bounds, fixed_vals, **kwargs
+    )
+
+
+def fixed_poi_fit(
+    poi_val, data, pdf, init_pars=None, par_bounds=None, fixed_params=None, **kwargs
+):
     r"""
     Run a maximum likelihood fit with the POI value fixed.
     This is done by minimizing the objective function of :func:`~pyhf.infer.mle.twice_nll`
@@ -142,6 +156,7 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs)
         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
+        fixed_params (`list`): Parameters to be held constant in the fit.
         kwargs: Keyword arguments passed through to the optimizer API
 
     Returns:
@@ -152,15 +167,11 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs)
         raise UnspecifiedPOI(
             'No POI is defined. A POI is required to fit with a fixed POI.'
         )
-    _, opt = get_backend()
-    init_pars = init_pars or pdf.config.suggested_init()
-    par_bounds = par_bounds or pdf.config.suggested_bounds()
-    return opt.minimize(
-        twice_nll,
-        data,
-        pdf,
-        init_pars,
-        par_bounds,
-        [(pdf.config.poi_index, poi_val)],
-        **kwargs,
-    )
+
+    init_pars = [*(init_pars or pdf.config.suggested_init())]
+    fixed_params = [*(fixed_params or pdf.config.suggested_fixed())]
+
+    init_pars[pdf.config.poi_index] = poi_val
+    fixed_params[pdf.config.poi_index] = True
+
+    return fit(data, pdf, init_pars, par_bounds, fixed_params, **kwargs)