Wrong poi_upper_limit computation with pyhf backend

[APMDSLHC]

System Settings

Fedora Linux 35 Python 3.9.12 spey 0.0.1

Describe the bug

The value of the poi upper limitreturned by poi_upper_limit() is not the one that gives a CLs of 0.05 when using the pyhf backend. This happens when get_multi_region_statistical_model() is called with a single SR.

To Reproduce

import spey, pyhf, jsonpatch, scipy

### Model definition
patch = [dict(
            op='add',
            path='/channels/0/samples/0',
            value=dict(
                name='sig',
                data=[0.4],
                modifiers=[
                    dict(
                        name='lumi',
                        type='lumi',
                        data=None
                    ),
                    dict(
                        name='mu_SIG',
                        type='normfactor',
                        data=None
                    )
                ]
            )
        )]

bkg = {
  "channels": [
    {
      "name": "SR1",
      "samples": [
        {
          "name": "bkg",
          "data": [
            0.8
          ],
          "modifiers": [
            {
              "data": None,
              "type": "lumi",
              "name": "lumi"
            }
          ]
        }
      ]
    }
  ],
  "measurements": [
    {
      "name": "BasicMeasurement",
      "config": {
        "poi": "mu_SIG",
        "parameters": [
          {
            "auxdata": [
              1
            ],
            "bounds": [
              [
                0.915,
                1.085
              ]
            ],
            "inits": [
              1
            ],
            "sigmas": [
              0.017
            ],
            "name": "lumi"
          }
        ]
      }
    }
  ],
  "observations": [
    {
      "name": "SR1",
      "data": [
        10
      ]
    }
  ],
  "version": "1.0.0"
}

### Getting the poi upper limit with spey
statModel = spey.get_multi_region_statistical_model("test",patch,bkg)
mu_ul_spey = statModel.poi_upper_limit(expected=spey.ExpectationType.observed,allow_negative_signal=True)
print(mu_ul_spey)

### Creation of the model
wsDict = jsonpatch.apply_patch(bkg, patch)
workspace = pyhf.Workspace(wsDict)
msettings = {'normsys': {'interpcode': 'code4'}, 'histosys': {'interpcode': 'code4p'}}
model = workspace.model(modifier_settings=msettings)

args={}
bounds = model.config.suggested_bounds()
bounds[model.config.poi_index] = [0,100]
args["par_bounds"] = bounds
args["return_expected"] = False
pver = float ( pyhf.__version__[:3] )
args["test_stat"] = "qtilde"

### Check if the upper limit on mu returned by spey is the correct one, i.e. is the one that gives a CLs of 0.05
CLs_spey = pyhf.infer.hypotest(mu_ul_spey, workspace.data(model), model, **args )
### It gives not 0.05 but 0.5
print(CLs_spey)

### Check if the problem really comes from poi_upper_limit
### Compute the poi upper limit via a simple script

### Define the function that must be 0 when the poi reaches its upper limit
def CLs_root_pyhf(mu, data, model, args):
    return 0.05 - float(pyhf.infer.hypotest(mu, data, model, **args))

data = workspace.data(model)
low = 0.1
high = 1.
while True:
    CLs_low = CLs_root_pyhf(low, data, model, args)
    CLs_high = CLs_root_pyhf(high, data, model, args)
    if CLs_low < 0 and CLs_high > 0:
        break
    elif CLs_low < 0 and CLs_high < 0:
        high *= 2.
    elif CLs_low > 0 and CLs_high > 0:
        low *= 1/2.
### Lower bound found for starting the Brent 's bracketing
print(CLs_low)
### Upper bound found for starting the Brent 's bracketing
print(CLs_high)
mu_ul_pyhf = scipy.optimize.brentq(lambda mu: CLs_root_pyhf(mu, data, model, args), low, high, rtol=1e-3, xtol=1e-3)
print(mu_ul_pyhf)

### Check if the upper limit on mu returned by pyhf backend only is the correct one, i.e. is the one that gives a CLs of 0.05
CLs_spey = pyhf.infer.hypotest(mu_ul_pyhf, workspace.data(model), model, **args )
### It gives 0.05 indeed
print(CLs_spey)

Expected behaviour

The poi upper limit computed by spey should be the one that gives a CLs of 0.05.

Additional information

The bounds for the Brent's bracketing are not the issue. The ones returned by spey (1.0, 64.0) are more or less the same as the ones found by the pyhf only script (0.1, 64.0). Replacing the former by the latter does not change the outcome.

[APMDSLHC]

Apparently it worked up to the commit 74a8f1c931dec1b02ea3495dbb84b3b89f094193 (Wed Feb 22 16:33:26).

[jackaraz]

Hi @APMDSLHC, I did some major changes in the code structure which fixed the issue on my end but please do all your checks again because these changes are big. So you can modify your code as follows

statModel = spey.get_multi_region_statistical_model("test", patch, bkg)
config = statModel.backend.model.config()
bounds = config.suggested_bounds
bounds[config.poi_index] = (config.minimum_poi, 100)
mu_ul_spey = statModel.poi_upper_limit(
    expected=spey.ExpectationType.observed, allow_negative_signal=False, par_bounds=bounds
)
print("muUL:", mu_ul_spey)
# Out: muUL: 38.38686452333036

Note that from now on, it won't do the adaptive bound adjusting; it seems like that is the thing that was giving you weird results. If you check the CLs distribution with respect to POI in the previous commit, the distribution has a discontinuity which caused by this adaptive bound rescaling. So I think we might need to come up with smarter ways to do the bound adjustment.

Please let me know if you find anything related to this.

[APMDSLHC]

Thanks for the fixes. Now I have an error when calling poi_upper_limit with SL backend:

File "/home/pascal/SModelS/spey/src/spey/optimizer/core.py", line 21, in fit
    print("model_configuration.fixed_poi_bounds():",model_configuration.fixed_poi_bounds())
AttributeError: 'list' object has no attribute 'fixed_poi_bounds'

The problem seems to arise because model_configuration is a list of zeros. Tell me if it's not a straightforward fix, I'll share you my code.

[jackaraz]

Thanks @APMDSLHC, I totally missed this one, I updated the necessary bits. It works on my side now:

stat_model_sl = spey.get_multi_region_statistical_model(
     "simple_sl_test",
     signal=[12.0, 11.0],
     observed=[51.0, 48.0],
     covariance=[[3.,0.5], [0.6,7.]],
     nb=[50.0, 52.0],
     delta_sys=0.,
     third_moment=[0.2, 0.1],
     xsection=0.5
 )
print(stat_model_sl.poi_upper_limit())
# Out: 0.8704245319083141

[APMDSLHC]

Thank you, it works well now, I am closing the issue.