Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: customizable confidence level for MINOS #497

Merged
merged 3 commits into from
Feb 2, 2025
Merged

feat: customizable confidence level for MINOS #497

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d23272c
allow changing confidence level for minos
alexander-held Feb 2, 2025
36dd421
add interface for fit functions
alexander-held Feb 2, 2025
35d666b
better comment
alexander-held Feb 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 23 additions & 3 deletions src/cabinetry/fit/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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(
Expand All @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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``.

Expand All @@ -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
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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.

Expand All @@ -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
Expand All @@ -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 = {}

Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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,
Expand Down
34 changes: 28 additions & 6 deletions tests/fit/test_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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] == {}


Expand Down Expand Up @@ -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
Expand All @@ -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] == {}


Expand All @@ -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,
Expand All @@ -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)],
Expand All @@ -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)],
Expand All @@ -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,
Expand All @@ -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)],
Expand All @@ -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)],
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand Down