Difference between uncorrelated and correlated shapes

[jackaraz]

Hi all, I have a very naive question; we are trying to understand how we should initialise our statistical model in different ways and need clarification. So let me set up an example;

Let's say I have a region with the following yields

signal_yields, bkg_yields, bkg_unc, obs_yields  = 1., 100., 1e-3, 100

Typically I would initialise this as

model = pyhf.simplemodels.uncorrelated_background([signal_yields], [bkg_yields], [bkg_unc])
data = [obs_yields] + model.config.auxdata

def computer(model, data):
    pars, twice_nllh = pyhf.infer.mle.fit(
        data,
        model,
        return_fitted_val=True,
        maxiter=200,
        par_bounds=model.config.suggested_bounds(),
    )
    return pars

And I would expect $\hat{\mu}$ to be zero since observed and bkg yields are the same; i.e.

print(computer(model, data))
# array([1.        , 1.00000002])

However, I got $\hat{\mu} = 1$. So I tested the same example in two ways; first, I used correlated background simple model

model = pyhf.simplemodels.correlated_background([signal_yields], [bkg_yields], [bkg_unc], [bkg_unc])
data = [obs_yields] + model.config.auxdata

pars = computer(model, data)
# pars = array([-7.03657102e-05,  1.73472348e-18])

Which is closer to what I expected than I tried adding 1 to bkg uncertainty in the uncorrelated simple model

model = pyhf.simplemodels.uncorrelated_background([signal_yields], [bkg_yields], [1+bkg_unc])
data = [obs_yields] + model.config.auxdata

pars = computer(model, data)
# pars = array([0.        , 0.99999957])

Which yields precisely what I expect. However, I don't know if I should add 1 to my uncertainty: is the uncertainty applied multiplicatively? Could you please confirm? Is this the same for correlated simple model? I guess for this particular example, since I have one bin, it doesn't matter, but I'm assuming that I should not be using a correlated simple model if I do not have upper and lower envelop information for the bkg uncertainty.

My second question is I see in the documentation that uncorrelated shape is always constrained via Poissonian. Is there a function that uses Gaussian instead for uncorrelated shapes? The reason I'm asking is that I'm trying to make the pdf similar to simplified likelihood.

Thanks!

[kratsg]

For the first part, I think before we answer a lot of the questions, you might want to refer to the documentation on the simplemodels:

• correlated_background
• uncorrelated_background

Already, I noticed that

model = pyhf.simplemodels.correlated_background([signal_yields], [bkg_yields], [bkg_unc], [bkg_unc])

is incorrect. What you'd want, assuming your bkg_unc = 1e-3 is an absolute uncertainty (and not relative? if it was relative, then bkg_unc = 0.1) is:

>>> model = pyhf.simplemodels.correlated_background([signal_yields], [bkg_yields], [bkg_yields - bkg_unc], [bkg_yields + bkg_unc])
>>> computer(model, data)
array([-0.00037948,  0.        ])

or if the background uncertainty was relative, then

>>> model = pyhf.simplemodels.correlated_background([signal_yields], [bkg_yields], [bkg_yields * (1 - bkg_unc)], [bkg_yields * (1 + bkg_unc)])
>>> computer(model, data)
array([3.45592019e-05, 2.41560547e-20])

since the correlated background is built by a histosys encoded with the up/down yields of the background sample. If you refer to the histosys documentation, there's an example that explains how this is meant to be encoded.

Likewise, for the uncorrelated_background, this is done with a shapesys (documentation) which is also not relative (e.g. not $\sigma_b / b$).

The other thing that's tripping you up here is that the order of the parameters is not consistent between these two models (and the auxdata changes too!)

>>> model.config.par_order
['mu', 'uncorr_bkguncrt']

and for correlated

>>> model.config.par_order
['correlated_bkg_uncertainty', 'mu']

so you need to probably (if easier) just change your computer to be

def computer(model, obs_yields):
    data = [obs_yields] + model.config.auxdata
    pars, twice_nllh = pyhf.infer.mle.fit(
        data,
        model,
        return_fitted_val=True,
        maxiter=200,
        par_bounds=model.config.suggested_bounds(),
    )
    return dict(zip(model.config.par_order, pars))

and making these changes, I get

>>> model = pyhf.simplemodels.uncorrelated_background([signal_yields], [bkg_yields], [bkg_unc])
>>> computer(model, obs_yields)
{'mu': 0.9999999999542389, 'uncorr_bkguncrt': 1.0000000207942246}

and

>>> model = pyhf.simplemodels.correlated_background([signal_yields], [bkg_yields], [bkg_yields - bkg_unc], [bkg_yields + bkg_unc])
>>> computer(model, obs_yields)
{'correlated_bkg_uncertainty': -0.00037948419484118976, 'mu': 0.0}

For the uncorrelated case, I suspect that the shapesys auxdata here is super larger that you end up with a weird edge case where the likelihood for $\mu = 1$ is slightly larger than the likelihood for $\mu = 0$:

>>> model.logpdf([0, 1], data)
array([-15.6542112])
>>> model.logpdf([1, 1], data)
array([-15.65917812])

but I have to think more about this here, but I think you might be setting your uncertainty incorrectly. @alexander-held ?

My second question is I see in the documentation that uncorrelated shape is always constrained via Poissonian. Is there a function that uses Gaussian instead for uncorrelated shapes? The reason I'm asking is that I'm trying to make the pdf similar to simplified likelihood.

See #820. Open issue atm. Constraint types cannot be changed from the config, but you can tweak the python code and do a bit of monkeypatching to make this work.

[alexander-held]

As far as I understand we are talking about this setup

signal_yields, bkg_yields, bkg_unc, obs_yields = 1.0, 100.0, 1e-3, 100
model = pyhf.simplemodels.uncorrelated_background(
    [signal_yields], [bkg_yields], [bkg_unc]
)
data = [obs_yields] + model.config.auxdata

which has the following spec:

{
    "channels": [
        {
            "name": "singlechannel",
            "samples": [
                {
                    "name": "signal",
                    "data": [1.0],
                    "modifiers": [{"name": "mu", "type": "normfactor", "data": null}],
                },
                {
                    "name": "background",
                    "data": [100.0],
                    "modifiers": [
                        {"name": "uncorr_bkguncrt", "type": "shapesys", "data": [0.001]}
                    ],
                },
            ],
        }
    ]
}

Single bin, 100 background events and 1 signal event at mu=1, tiny background uncertainty.

The NLLs from @kratsg look fine:

>>> model.logpdf([0, 1], data)
array([-15.6542112])
>>> model.logpdf([1, 1], data)
array([-15.65917812])

mu=1 has a smaller likelihood (set mu=100 or so to see it really drop down).

Of course the best-fit mu should be 0 for 100 observed events. I can reproduce the default wrong scipy result:

pyhf.infer.mle.fit(data, model)
# [1.         1.00000002]

Different init_pars do not help here, the model just refuses to fit anything. You can get a slightly more reasonable result with minuit:

pyhf.set_backend("numpy", "minuit")
pyhf.infer.mle.fit(data, model)
[0.05544482 1.        ]

One generic problem here is that we are at a parameter boundary for mu, you can relax that a bit with par_bounds=[(-10, 10), (-5, 5)] and that helps slightly. Beyond that you can set the tolerance to force further convergence:

pyhf.infer.mle.fit(data, model, par_bounds=[(-10, 10), (-5, 5)], tolerance=1e-5)
[-2.74080518e-06  9.99999998e-01]

There is something weird that happens with the small background uncertainty. Compare these two setups (again using scipy):

model = pyhf.simplemodels.uncorrelated_background([10], [100], [0.01])
data = [100] + model.config.auxdata
print(pyhf.infer.mle.fit(data, model))  # [0.99999993 0.99999993]

model = pyhf.simplemodels.uncorrelated_background([10], [100], [0.1])
data = [100] + model.config.auxdata
print(pyhf.infer.mle.fit(data, model))  # [4.4408921e-16 1.0000002e+00]

The only thing that changes is the background uncertainty. Likelihood scans look normal. Swapping out shapesys to staterror (thereby moving to a Gaussian constraint) also does not help, the spec for that is the following:

{
    "channels": [
        {
            "name": "singlechannel",
            "samples": [
                {
                    "name": "signal",
                    "data": [10],
                    "modifiers": [{"name": "mu", "type": "normfactor", "data": None}],
                },
                {
                    "name": "background",
                    "data": [100],
                    "modifiers": [
                        {"name": "uncorr_bkguncrt", "type": "staterror", "data": [0.01]}
                    ],
                },
            ],
        }
    ],
    "version": "1.0.0",
    "measurements": [
        {
            "name": "measurement",
            "config": {
                "poi": "mu",
                "parameters": [
                    {"bounds": [[0, 10]], "inits": [1.0], "fixed": False, "name": "mu"}
                ],
            },
        }
    ],
    "observations": [{"name": "singlechannel", "data": [100]}],
}

One last note about the original post: maxiter=200 may not work all that well generally, that is not a lot of steps to converge.

[jackaraz]

Thanks both; yes, I'm assuming my uncertainty is absolute, i.e. bkg_yields $ \pm $ bkg_unc. And should thank you both for the very detailed answer it clarified many things!

[alexander-held]

Correct, the number ending up in the shapesys "data" field is absolute.

[lukasheinrich]

Hi @jackaraz - I notice that the bkg uncertainty is very small - does it look better for e.g. a 5% uncertainty?

[alexander-held]

It works fine with larger uncertainties, but I have no idea what could break if they are small. Ideally we should recover the case with no background uncertainties just fine by making them smaller and smaller (minus likelihood offsets). I was wondering if we could run into any kind of under/overflow and that causes the issue?

[jackaraz]

Yes, as @alexander-held says, it works perfectly with larger unc values.