Decorrelate uncertainty for each bin

[Tomoya-Iizawa]

Hello,

I would like to ask if there is any way to implement systematics with decorrelated among bins?
For example, if we have bins [0, 250, 500, 750, 1000] and 2-point shape systematics as ttbar-ME, is it possible to decorrelate by each bin?

[alexander-held]

Hi, this is called a "shapesys" modifier in HistFactory and is supported by pyhf but not currently implemented in cabinetry as a config option. A part of the reason for that is that nobody had requested this modifier and its use is somewhat rare in my own experience. I created #454 to track this. The implementation should not be too difficult, there is currently just a somewhat long list of other things to also address (contributions here are certainly welcome).

The good news is that you can take an existing workspace and update it to decorrelate uncertainties in this manner. Here is an example snippet operating on a workspace called spec:

name_to_decorrelate = "Modeling"

for channel in spec["channels"]:
    for sample in channel["samples"]:
        if sample["name"] != "Background":
            continue  # veto samples as needed

        # find the normsys and histosys modifiers
        normsys = next(
            (
                mod
                for mod in sample["modifiers"]
                if mod["name"] == name_to_decorrelate and mod["type"] == "normsys"
            ),
            None,
        )
        histosys = next(
            (
                mod
                for mod in sample["modifiers"]
                if mod["name"] == name_to_decorrelate and mod["type"] == "histosys"
            ),
            None,
        )

        if normsys:
            norm_up, norm_down = normsys["data"]["hi"], normsys["data"]["lo"]
            sample["modifiers"].remove(normsys)  # delete the existing modifier
        else:
            norm_up, norm_down = 1.0, 1.0

        if histosys:
            # shapesys effect is symmetric, so symmetrize
            shapesys_up = np.asarray(histosys["data"]["hi_data"]) * norm_up
            shapesys_down = np.asarray(histosys["data"]["lo_data"]) * norm_down
            shapesys_data = np.abs(shapesys_up - shapesys_down) / 2
            sample["modifiers"].remove(histosys)  # delete the existing modifier
        else:
            shapesys_data = np.asarray(sample["data"]) * abs(norm_up - norm_down) / 2

        shapesys_modifier = {
            "data": shapesys_data.tolist(),
            "name": f"{name_to_decorrelate}_{channel['name']}_decorrelated_per_bin",
            "type": "shapesys",
        }
        sample["modifiers"].append(shapesys_modifier)

I have not tested this in a lot of detail but it seems to be doing the right thing. The idea is the following: define which modifiers you want to change (by name, veto samples as needed, can also veto channels if needed). Then find the modifiers that need to change: this assumes you have a NormPlusShape systematic in cabinetry, which by default would create a histosys modifier correlated with a normsys modifier. Then extract the required information from them, delete them from the modifier list, and lastly add a new modifier with the bin-by-bin decorrelated systematic uncertainty implemented as a shapesys.

[Tomoya-Iizawa]

Thank you for the answer and showing an example.

I have a few questions, one is why can this be bin-by-bin decorrelated? It looks like almost just removing normsys component.

Next question is that
shapesys_up = np.asarray(modifier["data"]["hi_data"]) * norm_up
shapesys_down = np.asarray(modifier["data"]["lo_data"]) * norm_down
Here I guess "modifier" should be "histosys", is it correct?

The last question is that even if I modified as above, it fails with error as
pyhf.exceptions.InvalidModel: Trying to add paramset shapesys/ttbar-ME-decorr on ttbar sample in HighMassOSs2thh_1bjets channel but other paramsets exist with the same name.
Do you have some idea how to resolve this error? I checked the modifiers list just before the line where the error is shown and ttbar-ME-decorr looks not duplicated in the HighMassOSs2thh_1bjets channel. This region is defined in second in the configuration next to HighMassOSs2thh_0bjets channel, so is it necessary to change the systematics source name for each channel (which sounds a little odd for me)? The log file is attached which shows the list of all the modeling uncertainties in the fit dumped just above the error line.
log2.txt

[alexander-held]

I have a few questions, one is why can this be bin-by-bin decorrelated? It looks like almost just removing normsys component.

It removes both histosys and normsys components (if they exist) to replace them by shapesys for the desired uncorrelated bin-by-bin treatment. The shapesys modifier itself has that behavior, it is controlled by independent parameters per bin (you'll notice more parameters in the fit results output).

shapesys_up = np.asarray(modifier["data"]["hi_data"]) * norm_up
shapesys_down = np.asarray(modifier["data"]["lo_data"]) * norm_down
Here I guess "modifier" should be "histosys", is it correct?

That was a bug, thanks for catching this! Indeed histosys is correct, I fixed it in the code above.

pyhf.exceptions.InvalidModel: Trying to add paramset shapesys/ttbar-ME-decorr on ttbar sample in HighMassOSs2thh_1bjets channel but other paramsets exist with the same name.

The reason is indeed that the name of the modifier needs to be unique per channel. I had overlooked that in the original implementation, fixed now also in the code above (by just making the channel name a part of the modifier name). The reason for this is that each bin in each channel is independent, and to keep them all independent the parameters per channel are identified with different names.

[Tomoya-Iizawa]

Thank you for the explanation and clarification

I tried implementing per channel as you indicated, but they are centered at 1 and very small errors, as the attached plot.

This is fit for asimov data, so I expected they should be centered on 0, and (hopefully) having +/- 1 sigma errors for systematics uncertainties.
Do you have some idea on this?

[alexander-held]

This is expected due to the different constraint term that is being used for shapesys parameters. See the pyhf documentation https://pyhf.readthedocs.io/en/v0.7.6/intro.html#id25, the relevant line in the table is "Uncorrelated Shape". The expected value for these parameters with Poisson constraint is 1 and the parameter uncertainty relates to the relative size of the impact of the parameter in each bin, which is usually a lot smaller than 1. In that sense they behave very differently in this plot when compared to the more common types of nuisance parameters with Gaussian constraints where the expectation is 0 +/- 1.

[Tomoya-Iizawa]

I see, thank you.
So for the case of the bin-by-bin uncorrelated uncertainty, my result looks as expected and the uncertainty becomes large if the size of impact of the parameter in each bin is large. But it seems to be difficult to compare the impact with other shape bin-by-bin correlated systematics, I guess comparing ranking would help to understand the impact.