pyhf: implementation of statistical uncertainty

user_14049393 0
Viewed 179

I have a question regarding the implementation of the statistical uncertainty. In the pyhf documentation https://scikit-hep.org/pyhf/likelihood.html#sample you mention that the way to infer statistical uncertainty is with modifier with "type": "staterror" and data field=[0.1].

So let's assume that I have a background channel which is coming from MC and I split my distribution in 3 bins:

"name": "background",
"data": [300., 50., 60.]

How do you correctly account for statistical uncertainty? From the construction of Poissonian pdf, I would say that by construction you already account for statistical uncertainty. Or do I have to have a modifier with staterror included? What exactly is the data field of staterror?

pyhf
1 Answers

I have a background channel which is coming from MC

Given that it sounds like you want to model the uncertainty in shape due to limited Monte Carlo sample size, the best modifier to use would be staterror. staterror is shared across all samples (that have a staterror modifier) in the bins it is applied with a Normal constraint, with the strength of the constraint being the per-sample uncertainties added in quadrature. Here the data key represents the absolute uncertainty in each bin of the sample (in this case being the Poisson uncertainty of the bin counts):

So, given your example of a single background sample with 3 bins, an example spec might look something like this (which I'll name bkg_only_spec.json)

{
    "channels": [
        {
            "name": "single_channel",
            "samples": [
                {
                    "name": "background",
                    "data": [
                        300.0,
                        50.0,
                        60.0
                    ],
                    "modifiers": [
                        {
                            "name": "uncorr_bkguncrt",
                            "type": "staterror",
                            "data": [
                                17.32051,
                                7.07107,
                                7.74597
                            ]
                        }
                    ]
                }
            ]
        }
    ],
    "observations": [
        {
            "name": "single_channel",
            "data": [
                300.0,
                50.0,
                60.0
            ]
        }
    ],
    "measurements": [
        {
            "name": "Measurement",
            "config": {
                "poi": "mu",
                "parameters": []
            }
        }
    ],
    "version": "1.0.0"
}

which we can see (note the constrained_by_normal) is still a valid spec with the CLI's inspect (though of course you need a signal sample as well to do any inference)

$ pyhf --version
pyhf, version 0.5.1
$ python answer.py
$ pyhf inspect bkg_only_spec.json
          Summary       
    ------------------  
       channels  1
        samples  1
     parameters  1
      modifiers  1

       channels  nbins
     ----------  -----
 single_channel    3  

        samples
     ----------
     background

     parameters  constraint              modifiers
     ----------  ----------              ----------
uncorr_bkguncrt  constrained_by_normal   staterror

    measurement           poi            parameters
     ----------        ----------        ----------
(*) Measurement            mu            (none)

where the below answer.py generates the spec.

# answer.py
import numpy as np
import json


def main():
    bins = [300.0, 50.0, 60.0]
    # rounding, as keeping full floating point is maybe a bit silly
    poisson_uncert = np.sqrt(bins).round(decimals=5).tolist()

    # just set the observations to be the same as the bin count here
    # as a placeholder
    spec = {
        "channels": [
            {
                "name": "single_channel",
                "samples": [
                    {
                        "name": "background",
                        "data": bins,
                        "modifiers": [
                            {
                                "name": "uncorr_bkguncrt",
                                "type": "staterror",
                                "data": poisson_uncert,
                            }
                        ],
                    }
                ],
            }
        ],
        "observations": [{"name": "single_channel", "data": bins}],
        "measurements": [
            {"name": "Measurement", "config": {"poi": "mu", "parameters": []}}
        ],
        "version": "1.0.0",
    }

    with open("bkg_only_spec.json", "w") as spec_file:
        json.dump(spec, spec_file, indent=4)


if __name__ == "__main__":
    main()
user_89319420
Related