diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index b4118da0..499bd6e9 100644 --- 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,6 +261,7 @@ 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, @@ -260,6 +269,7 @@ def _fit_model( 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,6 +342,7 @@ def _run_minos( minuit_obj: iminuit.Minuit, minos: Union[List[str], Tuple[str, ...]], labels: List[str], + minos_cl: Optional[float], ) -> Dict[str, Tuple[float, float]]: """Determines parameter uncertainties for a list of parameters with MINOS. @@ -337,6 +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) + minos_cl (Optional[float])): confidence level for the confidence interval Returns: Dict[str, Tuple[float, float]]: uncertainties indexed by parameter name @@ -347,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) + minuit_obj.minos(par_name, cl=minos_cl) minos_results = {} @@ -436,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, @@ -456,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, @@ -489,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 index c00d0651..93bb6e89 100644 --- 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,19 +322,25 @@ 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 level + 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() + # unknown parameter, MINOS does not run m = iminuit.Minuit(func_to_minimize, [1.0, 1.0]) 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", @@ -400,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, @@ -437,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, @@ -451,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, @@ -464,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,