Gaussian constraints on normfactor modifiers

[marfari]

In the current version of pyhf it is not possible to add Gaussian constraints on normfactor modifiers directly.

Let me know if I am wrong, but currently the only thing one can do is for e.g. given c = 10 +/- 2 to build a normfactor in the samples field:

{"data": None, "name": "c", "type": "normfactor"}

And to use a normsys modifier:

{"data": {"hi": 1+2/10, "lo": 1-2/10}, "name": "c_syst", "type": "normsys"}

And then to fix the normfactor in the parameters field:

{"name": "c", "bounds": [[0,1000]], "fixed": True, "inits": [10]}

The problem with this approach is when for e.g. c = 0 +/- 2 in which case the relative error becomes ill defined in the normsys modifier.

Is there a reason why a normfactor cannot have a Gaussian constraint term in the likelihood similar to what the lumi modifier has? Such that we could define "auxdata": [0] , "sigmas": [2] in the corresponding parameters field which is well defined in the case where the central value is 0?

In that case we would have in the samples field:

{"data": None, "name": "c", "type": "normfactor"}

And in the parameters field:

{"name": "c", "bounds": [[0, 1000]], "inits": [0], "auxdata": [0], "sigmas": [2]}

[alexander-held]

Hi @marfari, I think we might need a bit more information about the setup you have in mind. Once you fix the normfactor it no longer participates in your fits, so the first thing you describe to me sounds the same as using a normsys modifier and scaling your sample by a constant value of 10 in addition.

A normfactor is intended to be unconstrained by definition and the constrained version of it is a normsys modifier. The issue about negative predictions you point out is true, I would expect that this extends to lumi modifiers as well (I've not tested the behavior with them in practice for this case though).

One way to achieve something like a contribution of 0 +/- 2 events is to define a sample that provides a constant offset in each bin (e.g. by a value 3), assign a constant and negative normalization factor to it, and then re-define your sample which predicts 0 +/- 2 to instead predict e.g. 3 +/- 2 and when summed together with the constant negative samples your total prediction from those two parts becomes 0 +/- 2. You still need to ensure your total prediction from summing all other samples together does not become negative though since you otherwise will run into other issues evaluating the Poisson term for a negative predicted rate.

[marfari]

Hi @alexander-held ! Thank you for the reply!

Let me try to clarify a bit more my setup. I have 5 samples (signal, combinatorial background, and 3 physics background samples)

All my template histograms are normalised to 1 (so that the nominal rates are 1).

Then I have added the following modifiers (besides statistical error modifiers):
Signal: c*BF (normfactors) * c_syst (normsys)
Combinatorial: N_comb (normfactor)
3 phys. backgrounds: N_i (normfactor) * N_i_syst (normsys)

I have measurements for c and the three N_i: c +/- Delta_c and N_i +/- Delta_Ni; and I want to apply Gaussian constraints on these quantities.

The way I'm doing this right now is by having the normfactors c and N_i fixed (I'm scaling the nominal rates by these values) and then I have normsys with upward and downward fluctuations computed as 1+Delta_c/c, 1+Delta_Ni/N_i and 1-Delta_c/c, 1-Delta_Ni/N_i, respectively.

The issue is that sometimes the relative uncertainties can be larger or close to 100% or the N_i themselves can be 0 (e.g. 0+/- 2), which causes issues with the normsys modifier.

So I was just wondering wether it would make sense to be able to apply Gaussian constraints with the absolute uncertainties rather than with the relative uncertainties (similar to what the lumi modifier does) and if it could help with these corner cases.

I can try to do something like you suggest or to just remove these components from the fit when this happens, but just wondering if having this functionality could be helpful

[marfari]

To summarise: it would be nice if one could simply specify a normalisation modifier and be able to say I want to apply a Gaussian constraint on this normalisation which is centred at 0 and has an uncertainty of 2, using "auxdata": [0] and "sigmas": [2]

Currently this is not possible:

• normfactor does not support Gaussian constraints
• normsys assumes the nominal rate is 1 (so I need to multiple the shape by a normfactor which scales it by what I want and then fix this normfactor) and it fails when the central value is 0 because the relative error is not well defined. It also fails when the relative error is ~100% or larger.
• lumi does what I'd like but it can only be used once and it needs to be called "lumi", and I'd like to apply Gaussian constraints on several normalisations
• histosys would change the overall shape of the samples which is also not what I want

I'm not sure if the fit would complain if one applies a Gaussian constraint centred at 0 with uncertainty 2 (if one could use "auxdata": [0] and "sigmas": [2]) as you say

[alexander-held]

We probably need to distinguish conceptual physics points here from the question of how to implement something like you have in mind. Let's start with the physics: generally speaking it does not make sense for a sample to predict a negative amount of events in a bin (there can be some exceptions in practice e.g. for interference measurements depending on how they're set up). If you have a sample normalization estimated to be $0\pm2$, what does that mean? Generally I would read that as having an estimate consistent with 0 and an idea of how much larger than 0 it could be. What would make most sense to me for such a case is to implement an asymmetric uncertainty but enforce the prediction to not become negative. The normsys modifier does that through its use of exponential extrapolation, so you cannot use it around a central value of 0 (hence the suggestion to shift everything in my previous comment and then subtract it off again).

You can of course have a fit where you have a few samples contributing and you might have normfactor modifiers attached to them and then get a best-fit value for a normfactor of $0\pm2$. Importantly though in that fit there should be other samples contributing and keeping the total prediction per bin non-negative. The most elegant way to include that information in your fit would be to perform this auxiliary measurement in situ and have just a single fit that includes the setup you want to have plus an extra channel (or however many you might need) describing the likelihood that brought you to the estimate of $0\pm2$. That allows you to bypass any approximation and additionally it prevents you from neglecting any potential correlations between your estimated uncertainties that you would lose when manually propagating this from your N_i uncertainty estimates to your final fit.

Now on a purely technical level, here is what I think you would need to implement what you are describing: normsys will not work if you need your central value to be 0. You will need a histosys which linearly extrapolates. You do not have to encode a shape in your histosys, you can define it such that bin-by-bin the effect is the same. Take the background prediction you have, multiply it by the uncertainty of the corresponding N_i to define the histosys up and down variations, then set the nominal sample prediction to 0. Here is an example:

import pyhf

spec = {
    "channels": [
        {
            "name": "SR",
            "samples": [
                {
                    "data": [0],
                    "modifiers": [
                        {
                            "data": {"hi_data": [2], "lo_data": [-2]},
                            "name": "histosys_example",
                            "type": "histosys",
                        },
                    ],
                    "name": "Background",
                },
                {
                    "data": [10],
                    "modifiers": [
                        {
                            "data": None,
                            "name": "POI",
                            "type": "normfactor",
                        },
                    ],
                    "name": "Signal",
                }
            ],
        },
    ],
    "measurements": [{"config": {"parameters": [], "poi": "POI"}, "name": "example"}],
    "observations": [{"data": [10], "name": "SR"}],
    "version": "1.0.0",
}

ws = pyhf.Workspace(spec)
model = ws.model()
data = ws.data(model)

pyhf.set_backend("numpy", "minuit")
res = pyhf.infer.mle.fit(data, model, return_uncertainties=True)

for idx, r in enumerate(res):
    print(f"{model.config.par_order[idx]}: {r[0]:.3f} +/- {r[1]:.3f}")

While nothing spontaneously comes to mind that technically breaks this specific approach, I have never seen this idea used in practice personally. I would strongly suggest checking this further before using it, but really the better way to go in my opinion is to avoid this and find a different approach that does not run into the conceptual problems of predicting negative counts.

[marfari]

@alexander-held I see what you mean! I think for now I will remove the specific samples from the fit whenever the MC efficiency goes to zero. Thank you for the help!