add staterror config to spec to have faithful roundtrip

[lukasheinrich]

Description

we should add staterrorconfig as a channel property. Currently most analyses use Poisson as the config, which is not the default and thus needs to be specificed in the XML.

Converting XML -> JSON -> XML loses this information since we never record this non-default value in the JSON and thus cannot reproduce the original XML

during parsing, the staterrorconfig should become a value in the paramsets used by the staterror modifier.

[alexander-held]

@kratsg pointed out that ROOT 6.22 switches the HistFactory default constraint term to Poisson: https://root.cern/doc/v622/release-notes.html#release-6.2202

[alexander-held]

related: #662

[alexander-held]

Similarly to stat errors, ShapeSys has a ConstraintType option in HistFactory to switch to a Gaussian term. Presumably that information also is lost in a roundtrip.

[balunas]

I see this has been kicked back several versions, is there any estimate of when it might become available?

A related question: Do I understand correctly that if one wants to use a Poisson constraint, it can simply be specified as a shapesys rather than a staterror? As far as I can tell, the only difference aside from the Gaussian vs. Poisson is that staterror can be combined and pruned to reduce the total number of NPs.

[matthewfeickert]

I see this has been kicked back several versions, is there any estimate of when it might become available?

This is currently slated for v0.6.4. We're somewhat limited in our ability to give timelines for release schedules but if there is significant need for this to be implemented sooner the we can try to get a PR going once we have finished some currently outstanding work. We'd like to get this in within the next few weeks though for ourselves too. :slightly_smiling_face: The dev releases that are published to TestPyPI and master don't have official support or any promises of stability, but they are how we would recommend that people test releases in advance.

A related question: Do I understand correctly that if one wants to use a Poisson constraint, it can simply be specified as a shapesys rather than a staterror? As far as I can tell, the only difference aside from the Gaussian vs. Poisson is that staterror can be combined and pruned to reduce the total number of NPs.

While they are both multiplicative modifiers, in addition to the difference of the implementation between the Poisson and Normal constraint p.d.f. (c.f. the "Modifiers and Constraints" table in the docs, though it seems you already have), the method of determining the number of nuisance parameters that are allocated per sample for each modifier are different. So that's something to be aware of.

[balunas]

Thanks, that answers my question. It's not super urgent, since it's safe to use shapesys instead if I need a Poisson (the larger number of NPs is fine for my use case).

Out of curiosity, why is a Gaussian used for MC stat uncertainties? Since they're statistical in nature a Poisson should be more accurate, particularly in cases with large statistical uncertainties, right? The reason I got into asking these questions was because I was having issues with staterror in cases where the uncertainty is on the same scale as the nominal value. This only happens in a small number of bins, but it seems to give some questionable results where the Gaussian ~ Poisson approximation doesn't look accurate.

[balunas]

Just a gentle reminder on this, since it seems it's been pushed pack another major(?) version now. My analysis is now actually running into some limitations with the number of NPs, so it would be great to be able to combine Poisson staterror NPs via the Barlow-Beeston lite method which runs under the hood.

Again, why does staterror use a Gaussian constraint PDF at all? The Poisson is more correct and doesn't rely on large-stats assumptions. Is there any reason not to simply make that the default? I've seen several instances where it breaks down, forcing us to use shapesys, which then can't be combined into fewer NPs easily.

[alexander-held]

Hi @balunas, I can't comment on the time line, but maybe a few comments about the other points: I think the use of Gaussian as default came from the HistFactory implementation in ROOT, where this was the default choice before ROOT 6.22. I agree that the Poisson makes more sense as a default, at least when you have a uniform weight distribution.

If you do find that the choice of constraint term makes a big difference, then it may be worth revisiting the model setup. Large MC statistical uncertainties can bias your signal extraction (link to ATLAS-internal slides with a study). How are you building your workspace at the moment? I believe you should technically be able to correlate shapesys modifiers across all samples that you want to correlate to achieve what you are looking for. The correlation works by matching parameter names. This can also be done once you have a JSON workspace by renaming parameters as desired. I am probably overlooking some complication, but I currently do not see how shapesys makes this more complicated than staterror. There are some differences between the functionality of both approaches in ROOT (threshold-based pruning #662), but nothing I can spontaneously think of in the pyhf implementation.

If MC statistical uncertainties are problematic, it may also be interesting to consider re-binning or merging samples, which can address this type of issue more generally.

[balunas]

The drawback of just using shapesys is that, as far as I'm aware, the NPs for various contributions won't get combined into a single NP, resulting in unnecessarily slower fits. I'm aware that floating signal normalizations shouldn't be used while lumping signal MC stat uncertainties in with other NPs - this would of course result in incorrect scaling of the signal uncertainty contribution. But this isn't what I'm aiming to do. As long as you're careful about the signal treatment, there's no issue with working in a regime where the constraint function matters (and I wouldn't call these uncertainties problematic - in fact they don't always come from MC, but some other background estimation method which can have a statistical component from data). We've worked out a different solution for the time being, but for the future it would be really nice to be able to combine Poisson uncertainties this way.

There will probably be arguments against this, but I propose that Poisson be made the default for staterror in some future version.

[alexander-held]

the NPs for various contributions won't get combined into a single NP, resulting in unnecessarily slower fits

If you use the same parameter names, the modifiers are going to be correlated and you will end up with a single NP. I do not know whether there is internal optimization in pyhf that makes the interpolation/extrapolation calls faster with some custom staterror treatment though if you are referring to that. I'm not aware of anything like that, but maybe that exists.

Related to large MC uncertainties: this bias is not only an issue when scaling both signal and backgrounds together, but it also can appear with the more careful split treatment. When applying the same NP to both signal and background, the situation can certainly be worse though when pulls happen due to low background MC statistics and the effect is propagated to a (presumably high-statistic) signal.

I completely agree with you that Poisson should be the default, for the arguments you mention and for consistency with ROOT.

[matthewfeickert]

Just a gentle reminder on this, since it seems it's been pushed pack another major(?) version now.

Yes and no. We decided to not do a v0.6.4 patch release and just go for a v0.7.0 minor release (so it hasn't gotten pushed back, the release number changed). However, we're trying to get v0.7.0 done and out the door so that we can try to start having more frequent and useful patch and minor releases (we're at 5 months now since the last release but with a large amount of development that has happened since then — unshipped code is helping no one and I'm sure incentivizing people to just use untagged installs from GitHub which isn't good for reproducibility.)

I think the use of Gaussian as default came from the HistFactory implementation in ROOT, where this was the default choice before ROOT 6.22. I agree that the Poisson makes more sense as a default, at least when you have a uniform weight distribution.

What @alexander-held has here is correct. I don't think that there's any reasons to argue against changing it, so to help that along please make a new Issue for that.

[alexander-held]

Maybe to be practical, I think there are multiple aspects to this:

• change the default from Gaussian to Poisson,
• make the constraint term configurable,
• include the constraint term choice in the model specification.

Changing the default without implementing a way to switch behavior might not be ideal, as it will change the results of people who are used to the Gaussian setup. I think it might be possible to make this configurable though without already involving the JSON spec, and only do so in a separate step (assuming that this is maybe easier to make some partial progress). The natural way to configure this to me would be pyhf.Workspace.model via the modifier_settings kwarg, which already takes the interpolation codes for normsys and histosys, and could take another staterror key with something like {"constraint": "poisson"} or {"constraint": "gaussian"} as possible values.

[kratsg]

hi, i appreciate the discussion here. we do generally have limited time and I do think this is prioritized enough that we want to make sure it gets into the upcoming major release. What I propose we'll do (I have not discussed this with the core developers) is to support a user-friendly interface through pyhf in our next major release as this will change the API minimally, but allows for that support.

That said, it is possible to make this functional today in your existing code like so:

import pyhf

import json
ws = pyhf.Workspace(json.load(open('3b_tag21.2.27-1_RW_ExpSyst_79800_multibin_excl_Gtt_2400_5000_800.json')))

pdf = ws.model()
print(pdf.config.param_set('staterror_SR1L_Inj_Lmeff_cuts'))

def to_poisson(func):
    def wrapper(*args, **kwargs):
      result = func(*args, **kwargs)
      result['paramset_type'] = pyhf.parameters.constrained_by_poisson
      return result
    return wrapper

pyhf.modifiers.staterror.required_parset = to_poisson(pyhf.modifiers.staterror.required_parset)

pdf = ws.model()
print(pdf.config.param_set('staterror_SR1L_Inj_Lmeff_cuts'))

and this seems to work in being able to change things on a global-level (which I propose that pyhf supports primarily):

$ python do_poisson.py
<pyhf.parameters.paramsets.constrained_by_normal object at 0x1314b41f0>
<pyhf.parameters.paramsets.constrained_by_poisson object at 0x1317a0310>

This does also allow you to change constraints for other modifier types in a global way in a similar fashion ...

To summarize

• pyhf should support a global configuration for staterror constraint type (at least functionally in the code) which is an API change
• a future minor release can expand the v1.0.0 HiFa JSON specification to allow for this to be serializable
• the default staterror constraint should switch to Poisson to match ROOT 6.22+

[kratsg]

In the current HEAD of pyhf (pre-0.7.0), this works correctly

import pyhf
required_parset = pyhf.modifiers.staterror.required_parset

import json
ws = pyhf.Workspace(json.load(open('3b_tag21.2.27-1_RW_ExpSyst_79800_multibin_excl_Gtt_2400_5000_800.json')))

pdf = ws.model()
print(pdf.config.param_set('staterror_SR1L_Inj_Lmeff_cuts'))

def to_poisson(func):
    def wrapper(*args, **kwargs):
      result = required_parset(*args, **kwargs)
      result['paramset_type'] = 'constrained_by_poisson'
      result['factors'] = result.pop('sigmas')
      return result
    return wrapper

pyhf.modifiers.staterror.required_parset = to_poisson(pyhf.modifiers.staterror.required_parset)
pdf = ws.model()
print(pdf.config.param_set('staterror_SR1L_Inj_Lmeff_cuts'))

[alexander-held]

Related to roundtrips, an adversarial example with multiple staterror modifiers acting on the same sample & channel: #1830