Barlow-Beeston

[kskovpen]

Dear developers, could you please confirm whether pyhf supports the Barlow-Beeston implementation of statistical uncertainties? Apologies if I missed it in the documentation! Thanks a lot in advance, Kirill

[matthewfeickert]

@kskovpen (currently in v0.6.3) no. The HistFactory approach, as described in the original spec paper [CERN-OPEN-2012-016] (Section 2.2.1 Incorporating Monte Carlo statistical uncertainty on the histogram templates), is

Barlow and Beeston outlined a treatment of this situation in which each bin of each sample is given a nuisance parameter for the true rate, which is then fit using both the data measurement and the Monte Carlo estimate [?]. This approach would lead to several hundred nuisance parameters in the current analysis. Instead, the HistFactory employs a lighter weight version in which there is only one nuisance parameter per bin associated with the total Monte Carlo estimate and the total statistical uncertainty in that bin.

and then the rest of the section outlines things we mostly summarize in the docs.

Apologies if I missed it in the documentation.

Nothing to apologize about. While the MC Statistical Uncertainty (staterror) section of the modifiers docs does describe the implementation, as the section doesn't explicitly mention the phrase "Barlow-Beeston" it is easy to miss.

You might also want to look at Issue #760 for some general discussion.

[kskovpen]

Thanks @matthewfeickert ! One last follow-up question: was it ever considered to adjust the minimization procedure to use the analytic form of these parameters where one would just fix it to their best-fit values and remove it from the fit? Sorry again, if this is already available via some option :)

[matthewfeickert]

Probably, it's something along similar lines of #96

Correct.

[kskovpen]

Probably, it's something along similar lines of #96

[alexander-held]

Hi @kskovpen, it is possible to use full Barlow-Beeston in pyhf with some manual steps (depending on what you are starting with). The built-in staterror modifier results in BB "lite", and is intended to be used with all modifiers in a given channel having the same name¹.

We can use shapesys as a replacement for staterror, and assign one modifier per sample. There is one caveat with this: currently staterror has a Gaussian constraint (but see #1829), while shapesys has a Poisson constraint.

Here is a simple BB-lite example:

{
    "channels": [
        {
            "name": "SR",
            "samples": [
                {
                    "data": [25.0],
                    "modifiers": [
                        {
                            "data": [1.0],
                            "name": "staterror_SR",
                            "type": "staterror"
                        },
                        {"data": null, "name": "mu", "type": "normfactor"}
                    ],
                    "name": "Signal"
                },
                {
                    "data": [50.0],
                    "modifiers": [
                        {
                            "data": [3.0],
                            "name": "staterror_SR",
                            "type": "staterror"
                        }
                    ],
                    "name": "Background"
                }
            ]
        }
    ],
    "measurements": [
        {
            "config": {"parameters": [], "poi": "mu"},
            "name": "minimal_example"
        }
    ],
    "observations": [{"data": [75.0], "name": "SR"}],
    "version": "1.0.0"
}

Turning the staterror modifiers into shapesys just requires changing the type:

- "type": "staterror"
+ "type": "shapesys"

The two modifiers will still be correlated however, given that their names are the same.
To turn this into full Barlow-Beeston, you can make the modifier names per sample unique, for example like this:

"samples": [
    {
        "data": [25.0],
        "modifiers": [
            {
                "data": [1.0],
                "name": "staterror_SR_Signal",
                "type": "shapesys"
            },
            {"data": null, "name": "mu", "type": "normfactor"}
        ],
        "name": "Signal"
    },
    {
        "data": [50.0],
        "modifiers": [
            {
                "data": [3.0],
                "name": "staterror_SR_Background",
                "type": "staterror"
            }
        ],
        "name": "Background"
    }
]

Footnotes

1. Within pyhf that's technically not required, but will cause issues in ROOT and is probably best avoided (see https://github.com/scikit-hep/pyhf/issues/1830). ↩