Access modifier values

[ssaha1234]

Dear all,
I am writing a python script for pruning some of the systematics. I am able to access the name of the modifiers by using model.config.modifiers. I am now trying to figure out how to access the values of the modifiers (staterror, histosys, normsys) for pruning. Are there example scripts/syntax I can follow ?
Let me know if you need more details !

Thank you.
Best,
Shreya

[alexander-held]

Hi Shreya, what do you mean by "values of the modifiers"? Here is an example for how to access the constraint term width:

import pyhf

model = pyhf.simplemodels.uncorrelated_background([6], [9], [3])

for parameter in model.config.par_order:
    print(f"parameter {parameter} has type {dict(model.config.modifiers)[parameter]}")
    if model.config.param_set(parameter).constrained:
        print(f"  width of parameter: {model.config.param_set(parameter).width()}")

This might be useful for staterror pruning.

A nice (and I think quite general) approach may be using model._modifications, but I cannot comment much on that since I have not tried it yet. See #1652 and validate_systs.py.

A more convoluted method is getting nominal parameter settings for all parameters (e.g. cabinetry.model_utils.asimov_parameters), the associated pre-fit uncertainties (e.g. cabinetry.model_utils.prefit_uncertainties) and then for each parameter evaluating model.expected_data with the respective parameter varied up/down within uncertainties. That will return the yield for the full model, and taking the difference to the nominal yield could serve as a way to gauge the effect (see also model.main_model.expected_data with the return_by_sample keyword argument to get per-sample yields). This is most likely quite a bit less computationally efficient than model._modifications and involves a lot of duplicated calculations I think.

[ssaha1234]

Hi @alexander-held,
Thanks for the example and the detailed explanation ! I will read and apply it in my code and see how it goes. By the value of the modifiers, I meant something like this (for staterror, below) -

    "samples": [
        {
           "data": [26.0239, 20.4315, 3.28355], 
           "modifiers": [
              {"data": null, "name": "lumi", "type": "lumi"},
              {"data": [0.409128, 0.31449, 0.123071], "name": "TTZe_TTZeMCstat", "type": "staterror"},

I would like to access these values - [0.409128, 0.31449, 0.123071] for example, and similarly for histosys, normsys, etc.

Best,
Shreya

[ssaha1234]

Hi @alexander-held,
Apologies for the late reply ! I have written a very preliminary pruning code, specific to my analysis, I am happy to share it !

On a similar note, I have noticed that for staterror, where my Nominal value is 0, or staterror > 1 (which I set to 0, in my pruning code), I have an assertion error, when I try to access the model using pyhf. Is it something we expect ?

The error is as follows -
wspace = pyhf.Workspace(json.load(open("likelihood.json")))
model = wspace.model()

File "/afs/cern.ch/user/s/ssaha/.local/lib/python3.8/site-packages/pyhf/modifiers/staterror.py", line 116, in finalize assert len(sigmas[sigmas > 0]) == pdfconfig.param_set(mod).n_parameters AssertionError

Thanks !
Shreya

[alexander-held]

Hi @ssaha1234, the error you are seeing is (most likely) due to bins with zero predicted yield. See the discussion in #293. This now works on master as of #1657 and should be available in the next release as well. In addition, #1639 is probably interesting to link here as well, which will cause bins with zero staterror to be pruned automatically (also not included in a release yet). In the meantime, the only workaround I can spontaneously think of is assigning some very small non-zero prediction to these bins, or removing them from the model completely.

staterror > 1

Do you mean the relative size (staterror in model specification divided by sample yield), which can happen with negative weights present? If so, I am not sure why this should be pruned.

I have written a very preliminary pruning code, specific to my analysis, I am happy to share it !

I would be very interested in having a look! Feel free to share here, in a new post in this discussion forum, or in scikit-hep/cabinetry#311.

[ssaha1234]

Hi @alexander-held,
Yes, you are right ! I change the 0 Nominal yields and stat error values, to a very small number, and it seems to load fine, with no errors. I will use this as a workaround for now.

For the staterror > 1, I mean something like this, where the value of staterror is above > 1 (which I think is due to fluctuations and low yields in the region). As I understand, these values above 1 do not seem to have any physical meaning, please feel free to add your comments on it.
{"data": [1.50343, 1.28013, 0.834791], "name": "TTZe_TTZeMCstat", "type": "staterror"}

Sure, I can post the script in scikit-hep/cabinetry#311

Thanks a lot !
Shreya

[alexander-held]

The data fields for staterror modifiers contain the absolute statistical uncertainty, not the relative uncertainty. This admittedly is not clear from the documentation at the moment, see #1018. A 1% statistical uncertainty for a nominal yield of 1000 would result in a field with {"data": [10.0]}.

A relative uncertainty of >100% can also happen and be meaningful in the presence of negative MC weights. With large cancellations (due to negative weights), the absolute statistical uncertainty can be larger than the nominal yield.

[ssaha1234]

Yes, I assumed that the values here are the absolute statistical uncertainty and not the relative one. I am using hist.GetBinError(binIndex) to calculate it from my ntuples.

So following the example of 1000 events, the value of data for staterror can be greater than 1 and make sense. The numbers I am seeing in my likelihood file, I have Nominal data = 2.5 and staterror data is 1.5, so I have a 60% statistical uncertainty on this value.

For the relative uncertainty in terms of negative MC weights, my knowledge is limited, thanks for the explanation ! I have not seen absolute statistical uncertainty above the nominal yields yet, I will keep this in mind for future.

Thanks !
Best,
Shreya

[efemberg11]

@alexander-held Hi!
Im writing here, because I think my question is close for this topic.
When we adding or prunning our staterror modifiers, the amount of our NPs changes respectively. Can we prune one or couple statterror NP in model, like _prune_and_rename(.....)?
Also my question relate with TrexFitter's prune command 'MCstatThreshold', and I dont know how to realize it(
Correctly, how to control the number of gamma NPs in fit.

With regards, Ramdas

[alexander-held]

I don't think pyhf currently supports pruning a subset of bins, see #1143 for the related issue. Instead of pruning, you could hold the respective parameters (in the respective bins) constant in fits. See also #1279 which prevents you from being able to prune modifiers per sample / channel via the pyhf API. The setting provided by TRExFitter is used by HistFactory to drop nuisance parameters if the relative uncertainty is below the threshold, see #662. You could potentially achieve something comparable by parsing your pyhf model to determine which parameters to hold constant because they correspond to a term with an effect below this threshold (and then feed that list to any function performing fits).

[efemberg11]

Thank you for information!
Actually #662 is also my problem. As I understand from #662, there is no implementations, to prunning/adding gammas. And now, there is only manualy deleting/adding gamma_bins at spec "level"?
Also I mentioned, that when I reading spec with pyhf.readxml.parse(...), then making pyhf.Workspace(spec), then create model, data from workspace etc, and then I check my NP like this:
model.config.par_map['staterror_three_lep_presel_atLeast_3jets']
and i got
{'slice': slice(78, 91, None), 'paramset': <pyhf.parameters.paramsets.constrained_by_normal object at 0x7f260343ddf0>}

I have 13 gamma NPs (for channel 3_lep_atLeast_3jets), and they are constrained_by_normal
Does it mean, that when we reading workspace created by trexfitter with pyhf, we ignore whole this line:
StatErrorConfig RelErrorThreshold="0.05" ConstraintType="Poisson"
from our xml file?
With regards, Ramdas

[alexander-held]

And now, there is only manualy deleting/adding gamma_bins at spec "level"?

You can edit the spec, or prune full staterror parameters via the built-in pruning API / CLI in pyhf.

Does it mean, that when we reading workspace created by trexfitter with pyhf, we ignore whole this line
StatErrorConfig RelErrorThreshold="0.05" ConstraintType="Poisson"

Yes, this is correct and also leads to #760. The constraint term type is also ignored. You can (manually) replace staterror by shapesys to get Poisson constraints with current pyhf (0.7.0).