Computing p-value of an observed limit plot

[ssaha1234]

Dear pyhf users,
I am wondering if there is a provision to calculate the p-value of an observed limit, within pyhf ?
Thanks !
Best,
Shreya

[kratsg]

pyhf cls workspace.json reports the modified p-values $\mathrm{CL}_{\mathrm{S}}$ of the hypothesis test (for expected and observed). Is that what you're looking for? The python API is pyhf.infer.hypotest.

[matthewfeickert]

@ssaha1234 Can you be more specific about which $p$-value you want? If you've constructed an upper limit at the $(1-\alpha)$% confidence level then by construction you've found the value of your POI $\mu$ that results in the $\textrm{CL}_{s}$ (modified $p$-value) being $\alpha$. c.f. the pyhf.infer.hypotest and pyhf.infer.intervals.upperlimit docs for a super short bit on this as well as https://pyhf.github.io/pyhf-tutorial/HelloWorld.html#simple-upper-limit.

Did you want one of the $p$-values that appears in the constructing ratio of the $\textrm{CL}_{s}$ value?

$$
\mathrm{CL}{s} = \frac{\mathrm{CL}{s+b}}{\mathrm{CL}{b}} = \frac{p{s+b}}{1-p_{b}}
$$

If so, then use pyhf.infer.hypotest with return_tail_probs=True.

(All APIs in this answer given assuming using pyhf v0.6.3.)

[ssaha1234]

Hi @matthewfeickert and @kratsg,
Thanks for all the replies, and apologies for not giving much detail about the thing I was trying to achieve. Basically, I have run the hypotest over my json input file, testing 100 mu values (sys.argv[3] below) -

result = pyhf.infer.hypotest(float(sys.argv[3]), data, model, return_expected_set=True)

Now, I have already got the Brazil Band plot for all of my data (which shows the 95% expected limit on the branching fraction of the process vs different mass hypothesis of the new particle).
Next, I have made the observed line overlaying this CL expected Brazil Band using the same lines of code as above, and using unblinded data instead. As a next step, I would like to compute the observed signal significance as a function of the mass hypothesis. I do have the entire set of 5 values for the observed hypotest (-2, -1, 0, 1, 2)*sigma, and I am trying to get the significance from the result.
I think I need to return_tail_probs=True to get the p-values (and then the Z -values), let me know what you think !

Thanks a lot !
Best,
Shreya

[alexander-held]

Yes, maxiter has also been added to fit.significance, however there is no tagged version of cabinetry yet that will allow you to use it. The next version will be released when pyhf 0.7 comes out. You can however install the latest version directly from GitHub:

pip install --upgrade pyhf==0.7.0rc4
pip install --upgrade --force-reinstall git+https://github.com/scikit-hep/cabinetry.git@f0f9593

With that you can use maxiter in your cabinetry.fit.significance call.

Small point: you do not need to do all the imports in your script, import pyhf and import cabinetry is enough (unless you are using the other things you imported directly as well).

[alexander-held]

Alternatively you can call pyhf.infer.hypotest directly instead of using cabinetry and set the tolerance when setting the backend. See the cabinetry implementation for details. You can essentially copy that code into your script to make it work without having to upgrade any of your library versions.

[ssaha1234]

Hi @alexander-held,
Thanks for the points, I did something like this in my script now (added the hypotest instead of cabinetry) -

#!/usr/bin/env python3

import pyhf
from scipy.stats import norm
import json
import sys
import iminuit
import numpy as np

wspace = pyhf.Workspace(json.load(open('LH_pruned_unblinded_{}.json'.format(int(sys.argv[2])))))
model = wspace.model()
data = wspace.data(model)
 
init_pars = model.config.suggested_init()
par_bounds = model.config.suggested_bounds()
fixed_params = model.config.suggested_fixed()

obs_p_val, exp_p_val = pyhf.infer.hypotest(
    0.0,
    data,
    model,
    init_pars=init_pars,
    fixed_params=fixed_params,
    par_bounds=par_bounds,
    test_stat="q0",
    return_expected=True,
)
obs_p_val = float(obs_p_val)
exp_p_val = float(exp_p_val)
obs_significance = norm.isf(obs_p_val, 0, 1)
exp_significance = norm.isf(exp_p_val, 0, 1)

print('obs_p_val',obs_p_val)
print('exp_p_val',exp_p_val)
print('obs', obs_significance)
print('Exp', exp_significance)

I get the observed p-value as 0.5, which corresponds to an observed significance of 0. I am guessing, it means that something is wrong with the input file itself, or can there be some other issue, which I am not aware of (yet) ? I cannot find anything incorrect with the above lines of code.

Thanks !
Best,
Shreya

[alexander-held]

An observed significance of 0 happens when you do not see any signal. That doesn't necessarily mean that anything is wrong with your input. What is your best-fit signal strength? If it is <=0, then your significance is 0 (the calculation is capped and will not go below that for significance, corresponding to the p-value of 50%).

[ssaha1234]

Hi @alexander-held,
Thanks for the followup, I agree with it. I do not think the input is incorrect. I will check the value of the best-fit signal strength.
Thanks again for all the help !
Best,
Shreya