add interface for fit functions

scikit-hep · Feb 2, 2025 · 36dd421 · 36dd421
1 parent d23272c
commit 36dd421
Show file tree

Hide file tree

Showing 2 changed files with 46 additions and 14 deletions.
diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py
@@ -46,6 +46,7 @@ def _fit_model_pyhf(
     data: List[float],
     *,
     minos: Optional[Union[List[str], Tuple[str, ...]]] = None,
+    minos_cl: Optional[float] = None,
     init_pars: Optional[List[float]] = None,
     fix_pars: Optional[List[bool]] = None,
     par_bounds: Optional[List[Tuple[float, float]]] = None,
@@ -66,6 +67,8 @@ def _fit_model_pyhf(
         minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS
             algorithm for all parameters specified, defaults to None (does not run
             MINOS)
+        minos_cl (Optional[float]), optional): confidence level for the MINOS confidence
+            interval, defaults to None (use ``iminuit`` default of 68.27%)
         init_pars (Optional[List[float]], optional): list of initial parameter settings,
             defaults to None (use ``pyhf`` suggested inits)
         fix_pars (Optional[List[bool]], optional): list of booleans specifying which
@@ -117,7 +120,9 @@ def _fit_model_pyhf(
     best_twice_nll = float(best_twice_nll)  # convert 0-dim np.ndarray to float
 
     minos_results = (
-        _run_minos(result_obj.minuit, minos, labels) if minos is not None else {}
+        _run_minos(result_obj.minuit, minos, labels, minos_cl)
+        if minos is not None
+        else {}
     )
 
     fit_results = FitResults(
@@ -137,6 +142,7 @@ def _fit_model_custom(
     data: List[float],
     *,
     minos: Optional[Union[List[str], Tuple[str, ...]]] = None,
+    minos_cl: Optional[float] = None,
     init_pars: Optional[List[float]] = None,
     fix_pars: Optional[List[bool]] = None,
     par_bounds: Optional[List[Tuple[float, float]]] = None,
@@ -157,6 +163,8 @@ def _fit_model_custom(
         minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS
             algorithm for all parameters specified, defaults to None (does not run
             MINOS)
+        minos_cl (Optional[float]), optional): confidence level for the MINOS confidence
+            interval, defaults to None (use ``iminuit`` default of 68.27%)
         init_pars (Optional[List[float]], optional): list of initial parameter settings,
             defaults to None (use ``pyhf`` suggested inits)
         fix_pars (Optional[List[bool]], optional): list of booleans specifying which
@@ -234,7 +242,7 @@ def twice_nll_func(pars: np.ndarray) -> Any:
     corr_mat = m.covariance.correlation()  # iminuit.util.Matrix, subclass of np.ndarray
     best_twice_nll = m.fval
 
-    minos_results = _run_minos(m, minos, labels) if minos is not None else {}
+    minos_results = _run_minos(m, minos, labels, minos_cl) if minos is not None else {}
 
     fit_results = FitResults(
         bestfit,
@@ -253,13 +261,15 @@ def _fit_model(
     data: List[float],
     *,
     minos: Optional[Union[List[str], Tuple[str, ...]]] = None,
+    minos_cl: Optional[float] = None,
     init_pars: Optional[List[float]] = None,
     fix_pars: Optional[List[bool]] = None,
     par_bounds: Optional[List[Tuple[float, float]]] = None,
     strategy: Optional[Literal[0, 1, 2]] = None,
     maxiter: Optional[int] = None,
     tolerance: Optional[float] = None,
     custom_fit: bool = False,
+    cl: Optional[float] = None,
 ) -> FitResults:
     """Interface for maximum likelihood fits through ``pyhf.infer`` API or ``iminuit``.
 
@@ -274,6 +284,8 @@ def _fit_model(
         minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS
             algorithm for all parameters specified, defaults to None (does not run
             MINOS)
+        minos_cl (Optional[float]), optional): confidence level for the MINOS confidence
+            interval, defaults to None (use ``iminuit`` default of 68.27%)
         init_pars (Optional[List[float]], optional): list of initial parameter settings,
             defaults to None (use ``pyhf`` suggested inits)
         fix_pars (Optional[List[bool]], optional): list of booleans specifying which
@@ -300,6 +312,7 @@ def _fit_model(
             model,
             data,
             minos=minos,
+            minos_cl=minos_cl,
             init_pars=init_pars,
             fix_pars=fix_pars,
             par_bounds=par_bounds,
@@ -313,6 +326,7 @@ def _fit_model(
             model,
             data,
             minos=minos,
+            minos_cl=minos_cl,
             init_pars=init_pars,
             fix_pars=fix_pars,
             par_bounds=par_bounds,
@@ -328,8 +342,7 @@ def _run_minos(
     minuit_obj: iminuit.Minuit,
     minos: Union[List[str], Tuple[str, ...]],
     labels: List[str],
-    *,
-    cl: Optional[float] = None,
+    minos_cl: Optional[float],
 ) -> Dict[str, Tuple[float, float]]:
     """Determines parameter uncertainties for a list of parameters with MINOS.
 
@@ -339,8 +352,7 @@ def _run_minos(
         labels (List[str]]): names of all parameters known to ``iminuit``, these names
             are used in output (may be the same as the names under which ``iminiuit``
             knows parameters)
-        cl (Optional[float]), optional): confidence level for the confidence interval,
-            defaults to None (use ``iminuit`` default of 68.27%)
+        minos_cl (Optional[float])): confidence level for the confidence interval
 
     Returns:
         Dict[str, Tuple[float, float]]: uncertainties indexed by parameter name
@@ -351,7 +363,7 @@ def _run_minos(
             log.warning(f"parameter {par_name} not found in model")
             continue
         log.info(f"running MINOS for {par_name}")
-        minuit_obj.minos(par_name, cl=cl)
+        minuit_obj.minos(par_name, cl=minos_cl)
 
     minos_results = {}
 
@@ -440,6 +452,7 @@ def fit(
     data: List[float],
     *,
     minos: Optional[Union[str, List[str], Tuple[str, ...]]] = None,
+    minos_cl: Optional[float] = None,
     goodness_of_fit: bool = False,
     init_pars: Optional[List[float]] = None,
     fix_pars: Optional[List[bool]] = None,
@@ -460,6 +473,8 @@ def fit(
         minos (Optional[Union[str, List[str], Tuple[str, ...]]], optional): runs the
             MINOS algorithm for all parameters specified, defaults to None (does not run
             MINOS)
+        minos_cl (Optional[float]), optional): confidence level for the MINOS confidence
+            interval, defaults to None (use ``iminuit`` default of 68.27%)
         goodness_of_fit (bool, optional): calculate goodness of fit with a saturated
             model (perfectly fits data with shapefactors in all bins), defaults to False
         init_pars (Optional[List[float]], optional): list of initial parameter settings,
@@ -493,6 +508,7 @@ def fit(
         model,
         data,
         minos=minos,
+        minos_cl=minos_cl,
         init_pars=init_pars,
         fix_pars=fix_pars,
         par_bounds=par_bounds,

diff --git a/tests/fit/test_fit.py b/tests/fit/test_fit.py
@@ -101,7 +101,9 @@ def test__fit_model_pyhf(mock_minos, example_spec, example_spec_multibin):
 
     # including minos, one parameter is unknown
     model, data = model_utils.model_and_data(example_spec)
-    fit_results = fit._fit_model_pyhf(model, data, minos=["Signal strength", "abc"])
+    fit_results = fit._fit_model_pyhf(
+        model, data, minos=["Signal strength", "abc"], minos_cl=0.95
+    )
     assert fit_results.minos_uncertainty["Signal strength"] == (-0.1, 0.2)
     assert mock_minos.call_count == 1
     # first argument to minos call is the Minuit instance
@@ -110,6 +112,7 @@ def test__fit_model_pyhf(mock_minos, example_spec, example_spec_multibin):
         "Signal strength",
         "staterror_Signal-Region[0]",
     ]
+    assert mock_minos.call_args[0][3] == 0.95
     assert mock_minos.call_args[1] == {}
 
 
@@ -180,7 +183,9 @@ def test__fit_model_custom(mock_minos, example_spec, example_spec_multibin):
 
     # including minos
     model, data = model_utils.model_and_data(example_spec)
-    fit_results = fit._fit_model_custom(model, data, minos=["Signal strength"])
+    fit_results = fit._fit_model_custom(
+        model, data, minos=["Signal strength"], minos_cl=0.95
+    )
     assert fit_results.minos_uncertainty["Signal strength"] == (-0.1, 0.2)
     assert mock_minos.call_count == 1
     # first argument to minos call is the Minuit instance
@@ -189,6 +194,7 @@ def test__fit_model_custom(mock_minos, example_spec, example_spec_multibin):
         "Signal strength",
         "staterror_Signal-Region[0]",
     ]
+    assert mock_minos.call_args[0][3] == 0.95
     assert mock_minos.call_args[1] == {}
 
 
@@ -214,6 +220,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
     assert mock_pyhf.call_args[0][1] == data
     assert mock_pyhf.call_args[1] == {
         "minos": None,
+        "minos_cl": None,
         "init_pars": None,
         "fix_pars": None,
         "par_bounds": None,
@@ -228,6 +235,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
         model,
         data,
         minos=["Signal strength"],
+        minos_cl=0.95,
         init_pars=[1.5, 2.0],
         fix_pars=[False, True],
         par_bounds=[(0, 5), (0.1, 10.0)],
@@ -240,6 +248,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
     assert mock_pyhf.call_args[0][1] == data
     assert mock_pyhf.call_args[1] == {
         "minos": ["Signal strength"],
+        "minos_cl": 0.95,
         "init_pars": [1.5, 2.0],
         "fix_pars": [False, True],
         "par_bounds": [(0, 5), (0.1, 10.0)],
@@ -256,6 +265,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
     assert mock_custom.call_args[0][1] == data
     assert mock_custom.call_args[1] == {
         "minos": None,
+        "minos_cl": None,
         "init_pars": None,
         "fix_pars": None,
         "par_bounds": None,
@@ -270,6 +280,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
         model,
         data,
         minos=["Signal strength"],
+        minos_cl=0.95,
         init_pars=[1.5, 2.0],
         fix_pars=[False, True],
         par_bounds=[(0, 5), (0.1, 10.0)],
@@ -283,6 +294,7 @@ def test__fit_model(mock_pyhf, mock_custom, example_spec):
     assert mock_custom.call_args[0][1] == data
     assert mock_custom.call_args[1] == {
         "minos": ["Signal strength"],
+        "minos_cl": 0.95,
         "init_pars": [1.5, 2.0],
         "fix_pars": [False, True],
         "par_bounds": [(0, 5), (0.1, 10.0)],
@@ -310,15 +322,15 @@ def func_to_minimize(pars):
     m.errordef = 1
     m.migrad()
 
-    minos_results = fit._run_minos(m, ["b"], ["a", "b"])
+    minos_results = fit._run_minos(m, ["b"], ["a", "b"], minos_cl=None)
     assert np.allclose(minos_results["b"][0], -0.72622053)
     assert np.allclose(minos_results["b"][1], 0.47381869)
     assert "running MINOS for b" in [rec.message for rec in caplog.records]
     assert "b =  1.5909 -0.7262 +0.4738" in [rec.message for rec in caplog.records]
     caplog.clear()
 
     # custom confidence interval
-    minos_results = fit._run_minos(m, ["b"], ["a", "b"], cl=0.95)
+    minos_results = fit._run_minos(m, ["b"], ["a", "b"], minos_cl=0.95)
     assert np.allclose(minos_results["b"][0], -1.47037911)
     assert np.allclose(minos_results["b"][1], 0.81994008)
     caplog.clear()
@@ -328,7 +340,7 @@ def func_to_minimize(pars):
     m.errordef = 1
     m.migrad()
 
-    minos_results = fit._run_minos(m, ["x2"], ["a", "b"])
+    minos_results = fit._run_minos(m, ["x2"], ["a", "b"], minos_cl=None)
     assert minos_results == {}
     assert [rec.message for rec in caplog.records] == [
         "parameter x2 not found in model",
@@ -406,6 +418,7 @@ def test_fit(mock_fit, mock_print, mock_gof):
             (model, data),
             {
                 "minos": None,
+                "minos_cl": None,
                 "init_pars": None,
                 "fix_pars": None,
                 "par_bounds": None,
@@ -443,6 +456,7 @@ def test_fit(mock_fit, mock_print, mock_gof):
         (model, data),
         {
             "minos": None,
+            "minos_cl": None,
             "init_pars": init_pars,
             "fix_pars": fix_pars,
             "par_bounds": par_bounds,
@@ -457,10 +471,11 @@ def test_fit(mock_fit, mock_print, mock_gof):
     assert fit_results.bestfit == [1.0]
 
     # parameters for MINOS
-    fit_results = fit.fit(model, data, minos=["abc"])
+    fit_results = fit.fit(model, data, minos=["abc"], minos_cl=0.95)
     assert mock_fit.call_count == 3
     assert mock_fit.call_args[1] == {
         "minos": ["abc"],
+        "minos_cl": 0.95,
         "init_pars": None,
         "fix_pars": None,
         "par_bounds": None,
@@ -470,10 +485,11 @@ def test_fit(mock_fit, mock_print, mock_gof):
         "custom_fit": False,
     }
     assert fit_results.bestfit == [1.0]
-    fit_results = fit.fit(model, data, minos="abc", custom_fit=True)
+    fit_results = fit.fit(model, data, minos="abc", minos_cl=0.95, custom_fit=True)
     assert mock_fit.call_count == 4
     assert mock_fit.call_args[1] == {
         "minos": ["abc"],
+        "minos_cl": 0.95,
         "init_pars": None,
         "fix_pars": None,
         "par_bounds": None,