Extrapolation of fit results into validation regions

[miholzbo]

Hello,

Hopefully this is the correct way to ask such questions, otherwise please let my know in case there is a better communication channel!

I've recently started to play around with cabinetry, following the examples in cabinetry-tutorials; it's a really nice package! Some common application I couldn't find a way to do is to extrapolate the fit results to validation regions, i.e. regions not used in the fit. I think practically this means to just "copy over" the post-fit values (& uncertainties) of the normalization factors and nuisance parameters from a background-only fit of the control regions and evaluate the yields for such non-constraining regions (at least that's my understanding of what e.g. HistFitter does ..).

My naive approach was to define a second cabinetry configuration and call cabinetry.model_utils.prediction(model_VR, fit_results=fit_results) where fit_results comes from the fit of the CRs. But fit_results doesn't have the correct shape of course and also contains nuisance parameters related to the stat errors for the CRs which shouldn't be applied in the VRs.

Is there some way to achieve this, in particular to get post-fit plots of VRs? This could also be interesting to make post-fit plots of the CRs itself: often we define only 1-bin CRs for simplicity and then want to make post-fit plots of some kinematic variables in the CRs with several bins. Also here the "trick" in HistFitter was to define binned CRs as VRs to the fitted parameters would be evaluated there.

Any hint how to do this would be greatly appreciated (and apologies if I missed some description somewhere how to to this)!

Thanks!

Michael

[alexander-held]

Hi @miholzbo, this is a good place for such questions.

Your understanding of the basic approach is correct. As you note, the complication here is that the pyhf.Model which includes your validation regions might have additional parameters that are not part of your FitResults object (for example, MC statistical uncertainties for bins in the validation regions) and other parameters that are not needed for those validation regions. In addition, the parameter order might change, which also needs to be handled.

To help with this, there is a function model_utils.match_fit_results: you provide the target pyhf.Model (which contains the validation regions) and an existing FitResults object. It returns an updated version of that FitResults object that matches your target pyhf.Model, with reasonable assumptions for all the missing pieces of information that need to be filled in for that (see the documentation for what they are).

This conversion could in principle be done automatically within model_utils.prediction. The downside is that it hides away the assumptions going into this. I believe that it is useful to have users call model_utils.match_fit_results explicitly to make it very clear what is happening and give them more fine-grained control. I am always happy to hear opinions on the API design though and open to changes.

I opened an issue in the tutorial repository to add an example for this: cabinetry/cabinetry-tutorials#18. Given how common this use case is, I think this would be very useful.

cc @ekourlit (who has some ATLAS examples for this I believe)

[ekourlit]

Below an example workflow I am using for an ATLAS SUSY fit. Take care, this is for un-blinded scenarios. Otherwise, few extra steps should be taken. Hope this helps, in any case feel free to ask anything further.

# make a workspace containing three regions: SR1, VR1, CR1
ws = cabinetry.workspace.build(cabinetry_config)
# get model and observed data
model, data = cabinetry.model_utils.model_and_data(ws)
# remove SR1 and VR1 from the workspace
pruned_workspace = pyhf.Workspace(ws).prune(channels=['SR1', 'VR1'])
# get pruned model and pruned observed data
pruned_model, pruned_data = cabinetry.model_utils.model_and_data(pruned_workspace)
# fit 
fit_results_for_pruned_model = cabinetry.fit.fit(pruned_model, pruned_data)
# put back the SR1 and VR1 to the results by filling missing parameters with defaults
fit_results_for_full_model = cabinetry.model_utils.match_fit_results(model, fit_results_for_pruned_model)
# obtain post-fit model prediction
model_postfit = cabinetry.model_utils.prediction(model, fit_results=fit_results_for_full_model)

[miholzbo]

Hi @alexander-held,

Many thanks for the pointer to [model_utils.match_fit_results](https://cabinetry.readthedocs.io/en/latest/api.html#cabinetry.model_utils.match_fit_results), that works like a charm! I had looked at the API documentation of cabinetry but somehow must have browsed over this one .. I completely agree that this should be a dedicated step so people are aware of this subtlety.

@ekourlit, thanks a lot for sharing this snippet! The approach via pyhf.Workspace.prune() is indeed very neat and allows to work with only one cabinetry config.

Thanks again for the super-quick responses, my question has been perfectly addressed, so I'll close the discussion already.