shapefactor as parameter of interest

[TommyDESY]

Hi all,

I'm currently setting up a fit with 5 templates/samples. I included as a PoI a normfactor (Unconstrained Normalisation) which didn't cause any particular issue. I managed to create my pyhf model and managed to plot pulls, errors and fit with Asimov data.

However, what I actually want to do is use a shapefactor (Data-driven Shape). I use exactly the same code as previously but here an error arises. After building my pyhf spec I add the shapefactor as a modifier:

sample['modifiers'].append({'name': 'mu', 'type': 'shapefactor', 'data': None})

When I build my model with model = pyhf.Model(pyhf_spec), I get an AssertionError from pyhf/pdf.py l. 468:

467 s = self.par_slice(name)
--> 468 assert s.stop - s.start == 1

With a normfactor, this slice is indeed 1, but with shapefactor this slice has the same size as the number of bins in the fit (17 bins in lepton energy in my specific case).

I tried to remove this line from pdf.py but then this causes other errors when running the MLE fit. Therefore, it seems that this requirement is important.

What am I missing exactly here ? Is there anything specific that needs to be added to the pyhf spec when considering data-driven shapes ? Or is this due to some know issue ?

Thank you for your help.

[kratsg]

Hiya, so note that in pyhf, this limitation was at the moment due to the fact that we only allow for a single POI in a fit. When you have N POIs, it's not clear how you want to do the fit. See #179 for some of the discussion here. This is possibly a very good use-case for supporting shapefactor POIs (but this is a pretty uncommon one so far).

[alexander-held]

For the cases where the POI actually matters (limits / significance), I am not sure if you can make that work at the moment for an individual bin in that shapefactor modifier. At least I would not know how to achieve that using the existing pyhf API without re-implementing pieces. For everything else, setting the POI to None should be fine. The parameters won't technically be POIs, but the calculations are not affected by that definition. If you just want some data-driven estimate, then that sounds like you would not care about the POI label.

[TommyDESY]

Thank you, this is all very helpful. One more question on the technical side though.

When the POI is set to None, it doesn't seem possible to use minuit as an optimizer. In turn cabinetry cannot be used to perform the fit since minuit is automatically set as backend (despite the custom_fit option in cabinetry.fit.fit).

Again, what am I missing here ?

[alexander-held]

Could you please share (a simplified version of) the exact setup you are trying to use here? What you describe should generally work, but it is difficult to provide more specific advice without the details at this point.

[TommyDESY]

In my pyhf spec I have 1 channel (for now), 5 samples, 175 modifiers (nuisance parameters) + 1 shapefactor modifier. I want to fit 1 parameter in 17 bins (the shapefactor therefore corresponds to 17 modifiers in the spec).

I then build my model:
model = pyhf.Model(pyhf_spec, poi_name=None)

As backend I use jax and as optimizer minuit. However, this fails when I call the pyhf fit. I set custom_optimizer in set_backend to None for now. Then:

res = pyhf.infer.mle.fit(
    data=model.expected_data(model.config.suggested_init(),
                             include_auxdata=True), 
    pdf=model,
)

This works fine.

However, if I try with cabinetry it fails:

asimov_data = cabinetry.model_utils.asimov_data(model)
asimov_fit = cabinetry.fit.fit(model, asimov_data, custom_fit=True)

I get:
ValueError: Minimization failed, minimum is invalid.
from cabinetry/fit/__init__.py l.223

If I set custom_fit to False:
AttributeError: 'NoneType' object has no attribute 'correlation'
from pyhf/optimize/opt_minuit.py l. 134

I know that all this setup is fine if I simply use a normfactor instead of shapefactor as I've used it many times. Please, tell me if you want more information and thank you again for your help.

[alexander-held]

Please reduce this to a minimal setup you can share. The general functionality works just fine:

spec = {
    "channels": [
        {
            "name": "SR",
            "samples": [
                {
                    "data": [5.0, 10.0],
                    "modifiers": [
                        {"data": None, "name": "param", "type": "shapefactor"},
                    ],
                    "name": "signal",
                }
            ],
        }
    ],
    "measurements": [
        {
            "config": {"parameters": [], "poi": ""},
            "name": "test setup",
        }
    ],
    "observations": [{"data": [5, 10], "name": "SR"}],
    "version": "1.0.0",
}

import cabinetry
import pyhf

model, data = cabinetry.model_utils.model_and_data(spec)
cabinetry.fit.fit(model, data)

pyhf.infer.mle.fit(data, model)

The error you see is a minimization issue, which points to a problem in the model that you are using. That could come from many things, but I don't think it relates to the POI configuration.

[TommyDESY]

It turns out that the problem came from empty bins (obviously a problem when trying to rescale from 0). I avoided the issue by using fewer but larger bins.

In some cases I think it can also be fixed by adding a small positive epsilon to those empty bins but this might create other problems such as very large scaling factors (or hitting the boundary). In my case the minimisation still failed with this strategy.

I just write this here in case any one else stumbles upon this problem.

[matthewfeickert]

@TommyDESY This is now guarded against in pyhf v0.7.2 which was just released to PyPI.