How is pyhf dealing with asymmetric uncertainties ?

[TommyDESY]

Dear all,

In my current analysis (B --> Xu l nu at Belle II) we observe biases in the pull distributions obtained from toys. A minimal example of the pyhf_spec we're using is (1 channel, 3 samples, 1 POI, 1 NP):

{'channels': [{'name': 'Charged',
   'samples': [{'name': 'Xulnu',
     'data': [163.6619624355808,
      224.45580898062326,
      312.0552639523521,
      402.66760499449447,
      294.72046131570823,
      22.78226607432589],
     'modifiers': [{'name': 'mu_Xulnu', 'type': 'normfactor', 'data': None}]},
    {'name': 'Xclnu',
     'data': [2142.943246766925,
      2758.6685180701315,
      3097.5668878611177,
      2002.7195322178304,
      74.80800569057465,
      2.0612464249134064],
     'modifiers': [{'name': 'ff_DststBroad[0]',
       'type': 'histosys',
       'data': {'hi_data': [2184.1800931729376,
         2825.482837209478,
         3117.7314780633897,
         1996.440282246098,
         74.80800545215607,
         1.814018875360489],
        'lo_data': [2133.5764899738133,
         2743.8421267718077,
         3091.332954743877,
         2004.429974526167,
         74.80800545215607,
         2.1216768622398376]}},
      {'name': 'mu_Xclnu', 'type': 'normfactor', 'data': None}]},
    {'name': 'other_bkg',
     'data': [908.9859863962047,
      580.0524307717569,
      496.08645310113207,
      413.0712112658657,
      263.923080639448,
      44.730125417001545],
     'modifiers': [{'name': 'ff_DststBroad[0]',
       'type': 'histosys',
       'data': {'hi_data': [910.2731672744267,
         582.6203673486598,
         495.64627063134685,
         412.9712058152072,
         264.8207767768763,
         44.730125417001545],
        'lo_data': [908.6143180853687,
         579.4021521625109,
         496.13219153741375,
         413.09876992413774,
         263.7477066083811,
         44.730125417001545]}},
      {'name': 'mu_other_bkg', 'type': 'normfactor', 'data': None}]}]}]}

We produce a certain number of toys with:

model = pyhf.Model(pyhf_spec, poi_name='mu_Xulnu')
num_toys=5000
toys = model.make_pdf(
    pyhf.tensorlib.astensor(
        np.asarray(
            model.config.suggested_init(),
        )
    )
).sample((num_toys,))

We have 3 templates and for each we add a normfactor modifier. Then, the pulls are computed with:

param_index = model.config.par_names.index(param_name)
init_par = model.config.suggested_init()[param_index]
for toy in toys:
    fit_res = cabinetry.fit.fit(model, toy)
    mu_fit = fit_res.bestfit[param_index]
    uncert_fit = fit_res.uncertainty[param_index]
    pull = (mu_fit - init_par) / uncert_fit

When no nuisance parameters are added, the resulting pull distribution corresponds to a normal distribution as expected. However, we need to add histosys nuisance parameters to the fit and for some of those, the hi_data and lo_data yields are not symmetric around the nominal yield (both in their shape and total normalisation as in the spec above). You can see below that the Gaussian distribution is slightly shifted towards higher values of my POI.

My questions are therefore the following:
Is pyhf supposed to be able to deal with this kind of asymmetric errors ?
If yes, then what am I missing in the setup shown here ?
If no, then I suppose one needs to symmetrise these errors. We tried different ways which indeed eliminate the bias but force us to overestimate our systematic uncertainties. What would be the best way to do that in your opinion ?

Thank you in advance !

[alexander-held]

Hi @TommyDESY, pyhf can deal with asymmetric uncertainties just fine, but this can cause non-parabolic NLLs. You can have a look at the corresponding likelihood scan:

scan_results = cabinetry.fit.scan(model, data, "your_poi_name")
cabinetry.visualize.scan(scan_results)

which presumably is non-parabolic as well. I don't think there is anything wrong here and I do not see a need to symmetrize uncertainties. This asymmetric behavior can arise naturally (e.g. due to asymmetric systematic uncertainties or large systematic uncertainties acting on your signal).

[TommyDESY]

Thank you for the quick answer.

I see your point but then my real issue here is that I observe a bias in my signal extraction. Again, if my fit is working properly I would expect the pull distribution to be a normal distribution. Here, it seems that the asymmetric behaviour has an unexpected impact on the POI. Is there a way to mitigate that ? Or a different way to evaluate whether the fit is biased ?

[kratsg]

Typically, uncertainties with a high impact that are asymmetric will influence your fit that way. There's not much you can do to mitigate beyond symmetrizing and being more conservative with your uncertainties, but it's not (usually) an issue.