Fitting independent templates

[andrzejnovak]

Hi all, prompted by @matthewfeickert who I bothered originally :) I was wondering if I could use pyhf to find a best fit in setup where all templates are free floating. Here's what the problem could look like with code in "details". I suppose the place to start is the 2bin example, but the models seem to be built around the signal/background setup.

  import hist
  import numpy as np
  import mplhep as hep
  import matplotlib.pyplot as plt
  
  templates = []
  
  hall = hist.Hist(hist.axis.Regular(50, 0, 10))
  for i in [3, 5, 7]:
      h = hist.Hist(
              hist.axis.Regular(50, 0, 10),
          )
      h.fill(np.random.normal(i, 1, 10000))
      templates.append(h)
      
      hall.fill(np.random.normal(i, 1, round(10000*np.random.uniform(0.7, 1.3))))
      
  hep.set_style("CMS")
  hep.histplot(templates, stack=True, histtype='fill',label=["A", "B", "C"])
  hep.histplot(hall, histtype='errorbar', color='black', capsize=2, yerr=True, label='Total')
  plt.legend()

[alexander-held]

You can assign free-floating normfactor modifiers to all the three templates. This example can be converted into a workspace with one channel and three samples A, B, C. Each sample can have its own modifier, example for sample A:

"modifiers": [
    {
        "data": null,
        "name": "normalization_A",
        "type": "normfactor"
    }
]

Sample B would have a normfactor modifier with name "normalization_B" (or similar), similarly for C. If you only want to obtain the best-fitting model, it is irrelevant which of the normfactor modifiers you declare to be your parameter of interest in the model.

[alexander-held]

Here is a minimal example that you can add to your example code, it will fill in an appropriate workspace skeleton and perform a maximum likelihood fit:

import pyhf

spec = {
    "channels": [
        {
            "name": "some_region",
            "samples": [
                {
                    "data": [],
                    "modifiers": [
                        {"data": None, "name": "normalization_A", "type": "normfactor"}
                    ],
                    "name": "A",
                },
                {
                    "data": [],
                    "modifiers": [
                        {"data": None, "name": "normalization_B", "type": "normfactor"}
                    ],
                    "name": "B",
                },
                {
                    "data": [],
                    "modifiers": [
                        {"data": None, "name": "normalization_C", "type": "normfactor"}
                    ],
                    "name": "C",
                },
            ],
        }
    ],
    "measurements": [
        {
            "config": {"parameters": [], "poi": "normalization_A"},
            "name": "three_free_floating",
        }
    ],
    "observations": [{"data": [], "name": "some_region"}],
    "version": "1.0.0",
}

# fill in data
for i in range(3):
    spec["channels"][0]["samples"][i]["data"] = templates[i].values().tolist()
spec["observations"][0]["data"] = hall.values().tolist()

ws = pyhf.Workspace(spec)
model = ws.model()
data = ws.data(model)
fit_results = pyhf.infer.mle.fit(data, model)
for i, sample_name in zip(range(3), ["A", "B", "C"]):
    print(f"normfactor_{sample_name}: {fit_results[i]}")

[andrzejnovak]

Thanks a lot! It worked like a charm. Only hitch was I had to mod the template fetch with list(templates[i].to_numpy()[0])

If I understood the example correctly the scaling is given by "type": "normfactor", could this also take the form of an expression? Suppose this channel was split in two by some cut and the norm factor going up in one channel would have to be scaling the sample in the other channel down.

[alexander-held]

The .values() I was using might require a newer version of hist or boost-histogram (I used 2.3.0 and 1.0.1).

Expressions are not yet supported in pyhf, see #850.