Choosing alternative interpolator in normsys

[BlaiseDelaney]

Dear experts,

I wondered if there were a quick way to change the interpolator code in normsys. I understand that for releases >=6.0.0 the default interpolator changed. I am seeing some stark differences in a fit performed using pyhf v0.5.4 and v6.2.0 and I am relatively sure it's due to the normsys modifier. If possible, I would like to investigate whether the interpolator is the culprit.

In general, can I ask whether there is a clear criterion I should follow in choosing the best interpolator code? Apologies for my ignorance.

Thanks,

Blaise

[BlaiseDelaney]

Sorry, my bad. This was already answered pretty much here: in #1393.

Sorry for the noise!

[alexander-held]

Hi, you can set the interpolation code when creating the model:

workspace = pyhf.Workspace(spec)
model = workspace.model(
    modifier_settings={
        "normsys": {"interpcode": "code4"},
        "histosys": {"interpcode": "code4p"},
    }
)

The codes are described in the HistFactory manual, and the code4/code4p shown above are the defaults since pyhf 0.6. The available codes are listed in the Interpolators section in the API docs.

There is no "best" interpolation mechanism in general, but the defaults have some continuity in the derivatives that can be useful. The interpolation mechanism is ultimately a modeling choice, and ideally the dependence on this choice is small. If you see some very large differences, there may also be something else contributing.

[alexander-held]

cabinetry currently hardcodes these settings if you use model_utils.model_and_data, see the implementation. If switching out interpolation algorithms is commonly needed, we can add a flag for it (in my experience so far it is very rare to change from the default). You can however create your model and data directly with pyhf and continue using cabinetry with them. Have a look at the source code, it just uses pyhf directly to create the model and then there's a flag depending on the dataset you want.

[BlaiseDelaney]

Thank you for your reply. I wondered if I could also ask a related question: something I would really like to do is to translate the results I get from the fit for the normsys modifiers back to physical, constrained quantities. Can I ask if there is a quick way to do this? It would be very helpful to be able to think of the result as not some value constrained between ~[-1,1] but actually learn the best fit parameter and how it relates to the constraint I provide when creating the model specification.

Thanks a lot in advance.

[alexander-held]

Sorry, I do not understand what you are asking for. Maybe a good idea to open a new discussion since this topic sounds quite different from the topic here. That would allow others to find it more easily in the future.

To translate the result back into something more physically meaningful, you need to translate back from your simplified auxiliary measurement into whatever measurement you chose to approximate with it. That is very case dependent, and I do not see a general way to accomplish that. I might be misunderstanding the question though.

[alexander-held]

I think I might understand the question now: If you have a normsys modifier like

{
    "data": {
        "hi": 1.3,
        "lo": 0.8
    },
    "name": "par",
    "type": "normsys"
}

are you asking what value the associated sample is multiplied by if you pull par to some value? That is something for which pyhf most likely provides an API, but I do not know what it is by heart. One of the core devs might know, if you want to dig into the code I would suggest starting at pyhf.pdf.Model.expected_data which needs to calculate this somewhere further down to get the expected rates for given parameter values.

[BlaiseDelaney]

Thanks @alexander-held. What you suggest in your second reply is exactly what I have been doing so far. Apologies for the lack of clarity in my question.