diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index b4118da0..0029c4e5 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -367,7 +367,11 @@ def _run_minos( def _goodness_of_fit( - model: pyhf.pdf.Model, data: List[float], best_twice_nll: float + model: pyhf.pdf.Model, + data: List[float], + best_twice_nll: float, + *, + fix_pars: Optional[List[bool]] = None, ) -> float: """Calculates goodness-of-fit p-value with a saturated model. @@ -381,6 +385,8 @@ def _goodness_of_fit( data (List[float]): the observed data best_twice_nll (float): best-fit -2 log(likelihood) of fit for which goodness- of-fit should be calculated + fix_pars (Optional[List[bool]], optional): list of booleans specifying which + parameters are held constant, defaults to None (use ``pyhf`` suggestion) Returns: float: goodness-of-fit p-value @@ -416,7 +422,7 @@ def _goodness_of_fit( # of bins minus the number of unconstrained parameters n_dof = sum( model.config.channel_nbins.values() - ) - model_utils.unconstrained_parameter_count(model) + ) - model_utils.unconstrained_parameter_count(model, fix_pars=fix_pars) log.debug(f"number of degrees of freedom: {n_dof}") if n_dof <= 0: @@ -502,7 +508,9 @@ def fit( if goodness_of_fit: # calculate goodness-of-fit with saturated model - p_val = _goodness_of_fit(model, data, fit_results.best_twice_nll) + p_val = _goodness_of_fit( + model, data, fit_results.best_twice_nll, fix_pars=fix_pars + ) fit_results = fit_results._replace(goodness_of_fit=p_val) return fit_results diff --git a/src/cabinetry/model_utils.py b/src/cabinetry/model_utils.py index ec6533f8..a6a4cadc 100644 --- a/src/cabinetry/model_utils.py +++ b/src/cabinetry/model_utils.py @@ -481,8 +481,12 @@ def prediction( ) -def unconstrained_parameter_count(model: pyhf.pdf.Model) -> int: - """Returns the number of unconstrained parameters in a model. +def unconstrained_parameter_count( + model: pyhf.pdf.Model, + *, + fix_pars: Optional[List[bool]] = None, +) -> int: + """Returns the number of unconstrained, non-constant parameters in a model. The number is the sum of all independent parameters in a fit. A shapefactor that affects multiple bins enters the count once for each independent bin. Parameters @@ -490,15 +494,23 @@ def unconstrained_parameter_count(model: pyhf.pdf.Model) -> int: Args: model (pyhf.pdf.Model): model to count parameters for + fix_pars (Optional[List[bool]], optional): list of booleans specifying which + parameters are held constant, defaults to None (use ``pyhf`` suggestion) Returns: int: number of unconstrained parameters """ n_pars = 0 + # custom fixed parameters overrule suggested fixed parameters + if fix_pars is None: + fix_pars = model.config.suggested_fixed() + for parname in model.config.par_order: if not model.config.param_set(parname).constrained: # only consider non-constant parameters - n_pars += model.config.param_set(parname).suggested_fixed.count(False) + par_slice = model.config.par_slice(parname) + n_pars += fix_pars[par_slice].count(False) + return n_pars diff --git a/tests/fit/test_fit.py b/tests/fit/test_fit.py index c00d0651..16064201 100644 --- a/tests/fit/test_fit.py +++ b/tests/fit/test_fit.py @@ -353,7 +353,7 @@ def test__goodness_of_fit( assert mock_pars.call_args[1] == {} assert mock_count.call_count == 1 assert mock_count.call_args[0][0].spec == model.spec - assert mock_count.call_args[1] == {} + assert mock_count.call_args[1] == {"fix_pars": None} assert "Delta NLL = 0.084185" in [rec.message for rec in caplog.records] assert np.allclose(p_val, 0.91926079) caplog.clear() @@ -372,7 +372,7 @@ def test__goodness_of_fit( assert mock_pars.call_count == 2 # no new call (no auxdata) assert mock_count.call_count == 3 assert mock_count.call_args[0][0].spec == model.spec - assert mock_count.call_args[1] == {} + assert mock_count.call_args[1] == {"fix_pars": None} assert ( "cannot calculate p-value: 0 degrees of freedom and Delta NLL = 0.000000" in [rec.message for rec in caplog.records] @@ -480,7 +480,8 @@ def test_fit(mock_fit, mock_print, mock_gof): # goodness-of-fit test fit_results_gof = fit.fit(model, data, goodness_of_fit=True) - assert mock_gof.call_args_list == [((model, data, 2.0), {})] + assert mock_gof.call_args[0] == (model, data, 2.0) + assert mock_gof.call_args[1] == {"fix_pars": None} assert fit_results_gof.goodness_of_fit == 0.1 diff --git a/tests/test_model_utils.py b/tests/test_model_utils.py index b1369d6d..b74ef481 100644 --- a/tests/test_model_utils.py +++ b/tests/test_model_utils.py @@ -368,6 +368,12 @@ def test_unconstrained_parameter_count(example_spec, example_spec_shapefactor): model = pyhf.Workspace(example_spec_shapefactor).model() assert model_utils.unconstrained_parameter_count(model) == 3 + # fixed parameters can be provided to function + model = pyhf.Workspace(example_spec_shapefactor).model() + fix_pars = model.config.suggested_fixed() + fix_pars[0] = True + assert model_utils.unconstrained_parameter_count(model, fix_pars=fix_pars) == 2 + # fixed parameters are skipped in counting example_spec_shapefactor["measurements"][0]["config"]["parameters"].append( {"name": "Signal strength", "fixed": True}