Meaning of the expected CLs reported by pyhf.

[sabinekraml]

Dear pyhf developers,

We'd like to inquire about the exact definition of the expected CLs reported by pyhf.

As we understand it, expected limits, and thus also CLs_exp, are computed under the assumption that observation equals expectation (thus being agnostic of the actual observed data). We tried this out with the full likelihoods published by ATLAS, replacing in the json file the number of observed events in each signal region by the number of expected background events. This gives CLs_exp=CLs_obs, as it should be, but this value is different from the CLs_exp obtained with the original json file.

To give a concrete example from ATLAS-SUSY-2018-04 (hadronic stau search with two signal regions):

SR1cuts: #exp=6.057 #obs=10
SR2cuts: #exp=10.338 #obs=7

Here, #exp and #obs denote the number of expected SM events and the number of observed events, respectively.
Patching a signal of 1.568 events in region SR1cuts (=SRlow in the paper) and 7.701 events in region SR2cuts (=SRhigh in the paper), which corresponds to a benchmark point with m_stauLR = 400 GeV and m_LSP = 40 GeV, gives

CLs_obs = 0.94 and CLs_exp = 0.85.

However, when verifying the expected CLs "by hand", that is using

SR1cuts: #exp=6.057 #obs=6.057
SR2cuts: #exp=10.338 #obs=10.338

one obtains

CLs_obs = CLs_exp = 0.71 

for the same signal. Some clarification where the difference comes from would be much appreciated.

[lukasheinrich]

Dear Sabine The expected limit is the limit you would get if the observed data would be the yield you would expect for the background only hypothesis. Does that help? I can try giving an example Best Lukas

…

On Wed, 6 Oct 2021 at 17:59, Sabine Kraml ***@***.***> wrote: Dear pyhf developers, We'd like to inquire about the exact definition of the expected CLs reported by pyhf. As we understand it, expected limits, and thus also CLs_exp, are computed under the assumption that observation equals expectation (thus being agnostic of the actual observed data). We tried this out with the full likelihoods published by ATLAS, replacing in the json file the number of observed events in each signal region by the number of expected background events. This gives CLs_exp=CLs_obs, as it should be, but this value is different from the CLs_exp obtained with the original json file. To give a concrete example from ATLAS-SUSY-2018-04 (hadronic stau search with two signal regions): SR1cuts: #exp=6.057 #obs=10 SR2cuts: #exp=10.338 #obs=7 Here, #exp and #obs denote the number of expected SM events and the number of observed events, respectively. Patching a signal of 1.568 events in region SR1cuts (=SRlow in the paper) and 7.701 events in region SR2cuts (=SRhigh in the paper), which corresponds to a benchmark point with m_stauLR = 400 GeV and m_LSP = 40 GeV, gives CLs_obs = 0.94 and CLs_exp = 0.85. However, when verifying the expected CLs "by hand", that is using SR1cuts: #exp=6.057 #obs=6.057 SR2cuts: #exp=10.338 #obs=10.338 one obtains CLs_obs = CLs_exp = 0.71 for the same signal. Some clarification where the difference comes from would be much appreciated. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1619>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AARV6A2PSLUHGCK7KK4FKILUFRW7DANCNFSM5FO77KWQ> .

[sabinekraml]

Hi, this clarifies why using CLs_exp from pyhf does not give the result we expect, but it contradicts what we understand as expected values (used to evaluate, e.g., the sensitivity of an analysis). IMHO spelling this out in the documentation would be useful, as we (in this case the MadAnalysis5 team) are likely not the only ones who meet this pitfall.

[alexander-held]

The idea behind this approach is to make the "expected" result as realistic as possible, and that means e.g. taking into account possible pulls / non-nominal normalization factor values that would show up when confronted to real data. To be completely agnostic to real data, the "observed data" in the workspace has to be set to some pseudodata, and if set to Asimov then expected = observed limit by design.

[alexander-held]

Here is an example reproducing the "realistic" expected results manually. It is for significance instead of limit, but conceptually the same.

import pyhf

# Define a model where $\hat{\mu} \neq 1$ and pulls are present in best-fit:
model = pyhf.simplemodels.uncorrelated_background(
    signal=[12.0, 11.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
)
data = [70.0, 35.0] + model.config.auxdata

# MLE results:
def mle_fit(model, data):
    pyhf.set_backend("numpy", pyhf.optimize.minuit_optimizer(verbose=1))
    res = pyhf.infer.mle.fit(data, model, return_uncertainties=True)
    muhat, unc = res[model.config.poi_index]  # POI best-fit value and uncertainty
    print(f"mu = {muhat:.3f} +/- {unc:.3f}")

    theta_0, theta_unc_0 = res[1:][0]  # uncorr_bkguncrt, one per bin
    theta_1, theta_unc_1 = res[1:][1]
    print(f"NP 0: {theta_0:.3f} +/- {theta_unc_0:.3f}")
    print(f"NP 1: {theta_1:.3f} +/- {theta_unc_1:.3f}")


mle_fit(model, data)

# NP 1 is pulled -> "realistic" expected significance should take that into account.
# First calculate observed and expected significance via pyhf.infer.hypotest
obs_plus_exp = pyhf.infer.hypotest(
    0.0, data, model, test_stat="q0", return_expected=True
)
print(f"observed p-val: {obs_plus_exp[0]:.3%}, expected p-val {obs_plus_exp[1]:.3%}")

# To calculate expected significance, need to get auxdata from fit with $\mu=1$ fixed.
# Then build Asimov dataset with this auxdata -> MLE will have $\hat{\mu}=1$, and NPs
# will be at auxdata values.
fit_results = pyhf.infer.mle.fixed_poi_fit(1.0, data, model)
realistic_exp_data = model.expected_data(fit_results)

# Check MLE again. The best-fit result for $\mu$ is 1 by design. The NP values slightly
# differ from the unconstrained best-fit, since there $\hat{\mu}$ was different from 1.
# The pulls in this fit to the "realistic" Asimov dataset compensate that difference.
mle_fit(model, realistic_exp_data)

# Calculate realistic expected significance (by pretending the realistic Asimov dataset
# is observed data):
manual_realistic_exp = pyhf.infer.hypotest(
    0.0, realistic_exp_data, model, test_stat="q0"
)
print(f"expected p-val (manual calculation) {manual_realistic_exp:.3%}")

output:

mu = 0.407 +/- 0.504
NP 0: 1.043 +/- 0.061
NP 1: 0.810 +/- 0.096
observed p-val: 20.969%, expected p-val 2.942%
mu = 1.000 +/- 0.543
NP 0: 1.020 +/- 0.059
NP 1: 0.771 +/- 0.102
expected p-val (manual calculation) 2.942%

[lukasheinrich]

I'm curious @sabinekraml - what was the expectation (pun not intended) what the expected limits represent?

[matthewfeickert]

Here's @lukasheinrich's example from below as (overly formatted 😛) code:

# lukas_example.py
import pyhf

model = pyhf.simplemodels.uncorrelated_background(
    signal=[10], bkg=[50], bkg_uncertainty=[7]
)
observed_data = [55] + model.config.auxdata  # slight over fluctuation
expected_data_bkg_only = model.expected_data([0, 1])

init_params = model.config.suggested_init()
param_bounds = model.config.suggested_bounds()
fixed_params = model.config.suggested_fixed()

asimov_data_with_np = pyhf.infer.calculators.generate_asimov_data(
    0.0, observed_data, model, init_params, param_bounds, fixed_params
)

print(f"Background only model expected data: {expected_data_bkg_only}")
print(f"Asimov data:                         {asimov_data_with_np}")

test_poi = 1.0
# note the difference between "naive" bkg only and bkg only with data-derived pulls
print(f"\nCLs values for hypothesis of mu={test_poi}:\n")
print(f"{'Data set':29} | {'Observed CLs'} | {'Expected CLs'}")
print(f"{'-'*29} | {'-'*12} | {'-'*12}")
data_sets = [observed_data, expected_data_bkg_only, asimov_data_with_np]
data_descriptions = ["Observed data", "Background only expected data", "Asimov data"]
for data, description in zip(data_sets, data_descriptions):
    cls_observed, cls_expected = pyhf.infer.hypotest(
        test_poi, data, model, return_expected=True
    )
    print(f"{description:29} | {cls_observed:.7f}    | {cls_expected:.7f}")

$ python lukas_example.py 
Background only model expected data: [50.         51.02040816]
Asimov data:                         [52.47586877 53.54680487]

CLs values for hypothesis of mu=1.0:

Data set                      | Observed CLs | Expected CLs
----------------------------- | ------------ | ------------
Observed data                 | 0.4541889    | 0.3279657
Background only expected data | 0.3163970    | 0.3163970
Asimov data                   | 0.3279684    | 0.3279708

[kratsg]

Just a quick comment for now -- but the computation of CLs requires first - a fit to generate Asimov data - and this depends on the observed data you specify. This Asimov dataset will depend on the observed dataset. Does that clarify?

[kratsg]

One thing @alexander-held reminded me of to clarify, that goes along with what @lukasheinrich said. The Asimov will center your background expectation to the observed data.

If you switch from observed data to expected data -- your background-only model will have NPs centered at different values because they're fit to different datasets. The expected CLs has to be different.

[lukasheinrich]

Here is an example: [image: image.png]

…

On Wed, Oct 6, 2021 at 6:09 PM Giordon Stark ***@***.***> wrote: One thing @alexander-held <https://github.com/alexander-held> reminded me of to clarify, that goes along with what @lukasheinrich <https://github.com/lukasheinrich> said. The Asimov will center your background expectation to the observed data. If you switch from observed data to expected data -- your background-only model will have NPs centered at different values because they're fit to different datasets. The expected CLs has to be different. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1619 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AARV6A23QS3ZRWKQXS5ICD3UFRYDFANCNFSM5FO77KWQ> .

[sabinekraml]

the image doesn't show ....

[matthewfeickert]

Attachments from email responses don't get rendered. :/ When Lukas can get back to the GitHub web UI he can post the example image in the thread of his original response above and we can delete this comment section.

[lukasheinrich]

a apologies - here. .note the differenec between "naive" b-only and b-only with data-derivd pulls

[sabinekraml]

Hi guys, thanks for the detailed answers. So pyhf gives an "a-posteriori" expected CLs, and this can be quite different from the "a-priori" one where one sets data=SM (e.g. before unblinding and/or to determine sensitivities). I think it is worthwhile to make this very clear for the user not only by quoting the underlying formulae, like in #1367, but giving also some physics understanding. The point is that many of the official limit plots from the experiments show the a-priori expected limit, not the a-posteriori one. So when using a published full likelihood from ATLAS and naively comparing pyhf's CLs_exp with the expected limit from the analysis, there is a discrepancy. Perhaps experimentalists are all aware of this, but it is a potential serious pitfall for users from the theory side (recasters.....).

[matthewfeickert]

I think it is worthwhile to make this very clear for the user not only by quoting the underlying formulae, like in #1367, but giving also some physics understanding.

We're quite happy to try and make this more clear for people. So we'll address this in:

• Issue Clarify in docstrings what an expected CLs value is #1620
• Add example of what the expected and observed CLs values actually are pyhf/pyhf-tutorial#10

The point is that many if not most of the official limit plots from the experiments show the a-priori expected limit, not the a-posteriori one. So when using a published full likelihood from ATLAS and naively comparing pyhf's CLs_exp with the expected limit from the analysis, there is a discrepancy.

Can you point us to some examples of this @sabinekraml? I think @lukasheinrich, @kratsg, @alexander-held, and I are all surprised to hear this as the ATLAS SUSY group (and I'm pretty sure the ATLAS Higgs groups as well) definitely uses Asimov data to produce expected limits (what you refer to as "a-posteriori"). This is also what is recommended in the joint recommendations from ATLAS and CMS in Procedure for the LHC Higgs boson search combination in Summer 2011 Section 2 (there the phrase "pseudo-data" is used) and (using "Asimov") in the February, 2011 ATLAS Frequentist Limit Recommendation Section 3.

[edit: I had misremembered what the Procedure for the LHC Higgs boson search combination in Summer 2011 actually recommends (sorry for not reading before posting) and it isn't giving a prescription that uses something like the Asimov dateset. So please just refer to the 2011 ATLAS Frequentist Limit Recommendation for this question.]

Perhaps experimentalists are all aware of this, but it is a potential serious pitfall for users from the theory side (recasters.....).

Yes, it is clear that there is not good inter-field communication on what is happening here. This is at least good to know about so that we can try to avoid and fix this in the future. Thank you for raising this though, as it is quite important!

[sabinekraml]

OK, maybe I'm mistaken here; will get back to the books. (But, pseudo-data can also be sampled from the SM hypothesis without knowledge of the measured data, can't it) In any case, my original question was not about the difference between expected and observed limits as on pyhf/pyhf-tutorial#10, but on the different ways to compute an expected limit, and how this is done in pyhf. This is clarified. What we have to do on our side is to use one method consistently when recasing, and be doubly careful when comparing with expected limits in the experimental publications, in order to compare apples with apples.

…

On Thu, Oct 7, 2021 at 11:13 PM Matthew Feickert ***@***.***> wrote: Hi guys, thanks for the detailed answers. So pyhf gives an "a-posteriori" expected CLs, and this can be quite different from the "a-priori" one where one sets data=SM (e.g. before unblinding and/or to determine sensitivities). I think it is worthwhile to make this very clear for the user not only by quoting the underlying formulae, like in #1367 <#1367>, but giving also some physics understanding. The point is that many if not most of the official limit plots from the experiments show the a-priori expected limit, not the a-posteriori one. So when using a published full likelihood from ATLAS and naively comparing pyhf's CLs_exp with the expected limit from the analysis, there is a discrepancy. Perhaps experimentalists are all aware of this, but it is a potential serious pitfall for users from the theory side (recasters.....). I think it is worthwhile to make this very clear for the user not only by quoting the underlying formulae, like in #1367 <#1367>, but giving also some physics understanding. We're quite happy to try and make this more clear for people. So we'll address this in: - Issue #1620 <#1620> - pyhf/pyhf-tutorial#10 <pyhf/pyhf-tutorial#10> The point is that many if not most of the official limit plots from the experiments show the a-priori expected limit, not the a-posteriori one. So when using a published full likelihood from ATLAS and naively comparing pyhf's CLs_exp with the expected limit from the analysis, there is a discrepancy. Can you point us to some examples of this @sabinekraml <https://github.com/sabinekraml>? I think @lukasheinrich <https://github.com/lukasheinrich>, @kratsg <https://github.com/kratsg>, @alexander-held <https://github.com/alexander-held>, and I are all surprised to hear this as the ATLAS SUSY group definitely uses Asimov data to produce expected limits. This is also what is recommended in the joint recommendations from ATLAS and CMS in Procedure for the LHC Higgs boson search combination in Summer 2011 <https://inspirehep.net/literature/1196797> Section 1 (there the phrase "pseudo-data" is used) and (using "Asimov") in the February, 2011 ATLAS Frequentist Limit Recommendation <https://indico.cern.ch/event/126652/contributions/1343592/attachments/80222/115004/Frequentist_Limit_Recommendation.pdf> Section 3. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1619 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG3ROOG277UFIJDYVCVYF5LUFYENXANCNFSM5FO77KWQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

--

_______________________________________________________ Sabine Kraml - ***@***.*** - www.kraml.net LPSC 53 Av des Martyrs 38026 Grenoble France (+33)(0)4 76 28 40 52

[matthewfeickert]

(But, pseudo-data can also be sampled from the SM hypothesis without knowledge of the measured data, can't it)

Yeah, I misremembered what was in that recommendation and didn't' re-read before posting. I've edited my #1619 (reply in thread) above.

In any case, my original question was not about the difference between expected and observed limits as on pyhf/pyhf-tutorial#10, but on the different ways to compute an expected limit, and how this is done in pyhf. This is clarified.

Right. I want to address both of those things (how both are calculated and the prescription used for expected limits in ATLAS) in the that Issue. 👍

What we have to do on our side is to use one method consistently when recasting, and be doubly careful when comparing with expected limits in the experimental publications, in order to compare apples with apples.

I don't actually know what CMS does in these cases given my edited response above, so that will be worth finding out as well.

[matthewfeickert]

I don't actually know what CMS does in these cases given my edited response above, so that will be worth finding out as well.

The good news is that CMS Combine is indeed doing the same thing that ATLAS does (e.g. using the Asimov data for expected limits) when they are performing Frequentist limits (they may also set Bayesian limits, but that is apparently not common these days).

In the Combine docs, there is a section on Asymptotic Frequentist Limits that confirms that they use the Asimov approach as well

The AsymptoticLimits calculation follows the frequentist paradigm for calculating expected limits. This means that the routine will first fit the observed data, conditionally for a fixed value of r and set the nuisance parameters to the values obtained in the fit for generating the Asimov data, i.e it calculates the post-fit or a-posteriori expected limit.

This is doubly confirmed when looking at the Combine FAQ "Why does changing the observation in data affect my expected limit?"

The expected limit (if using either the default behaviour of -M AsymptoticLimits or using the LHC-limits style limit setting with toys) uses the post-fit expectation of the background model to generate toys. This means that first the model is fit to the observed data before toy generation.

Though it is worth noting that Combine has a "Blind limits" option

In order to use the pre-fit nuisance parameters (to calculate an a-priori limit), you must add the option --noFitAsimov or --bypassFrequentistFit.

Though I would assume (read as: someone who knows should correct me) that CMS would make this very explicit if they were to ever present limits like this given the Asimov procedure being viewed as the "default" in ATLAS and CMS.

(Thanks @bregnery for pointing me in the right direction on the Combine docs. 👍)

[sabinekraml]

In the Combine docs, there is a section on Asymptotic Frequentist Limits

<https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part3/commonstatsmethods/#asymptotic-frequentist-limits> that confirms that they use the Asimov approach as well Yes, and scrolling down to "Computing the expected significance" one finds the "a-priori expected" as the default behavior. Anyways, it would be really good to discuss "live" if you have time next week.

…

On Fri, Oct 8, 2021 at 4:19 PM Matthew Feickert ***@***.***> wrote: I don't actually know what CMS does in these cases given my edited response above, so that will be worth finding out as well. The good news is that CMS Combine is indeed doing the same thing that ATLAS does (e.g. using the Asimov data for expected limits) *when* they are performing Frequentist limits (they may also set Bayesian limits, but that is apparently not common these days). In the Combine docs, there is a section on Asymptotic Frequentist Limits <https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part3/commonstatsmethods/#asymptotic-frequentist-limits> that confirms that they use the Asimov approach as well The AsymptoticLimits calculation follows the frequentist paradigm for calculating expected limits. This means that the routine will first fit the observed data, conditionally for a fixed value of r and set the nuisance parameters to the values obtained in the fit for generating the Asimov data, i.e it calculates the *post-fit* or *a-posteriori* expected limit. This is doubly confirmed when looking at the Combine FAQ <http://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part4/usefullinks/#faq> "Why does changing the observation in data affect my expected limit?" The expected limit (if using either the default behaviour of -M AsymptoticLimits or using the LHC-limits style limit setting with toys) uses the *post-fit* expectation of the background model to generate toys. This means that first the model is fit to the *observed data* before toy generation. Though it is worth noting that Combine has a "Blind limits" option <https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part3/commonstatsmethods/#blind-limits> In order to use the *pre-fit* nuisance parameters (to calculate an *a-priori* limit), you must add the option --noFitAsimov or --bypassFrequentistFit. Though I would assume (read as: someone who knows should correct me) that CMS would make this *very* explicit if they were to ever present limits like this given the Asimov procedure being viewed as the "default" in ATLAS and CMS. (Thanks @bregnery <https://github.com/bregnery> for pointing me in the right direction on the Combine docs. 👍) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1619 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AG3ROODKQJ3G2KIC4VG7N2LUF34YJANCNFSM5FO77KWQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

--

_______________________________________________________ Sabine Kraml - ***@***.*** - www.kraml.net LPSC 53 Av des Martyrs 38026 Grenoble France (+33)(0)4 76 28 40 52

[kpedro88]

In my experience in CMS, limits or significances from published analyses always use the Asimov approach to get the expected results. The blind/a-priori behavior is useful for optimizing a blinded analysis before unblinding, as indicated at https://cms-analysis.github.io/HiggsAnalysis-CombinedLimit/part3/commonstatsmethods/#computing-the-expected-significance (emphasis added):

a-priori expected (the default behavior): does not depend on the observed dataset, and so is a good metric for optimizing an analysis when still blinded.

[WolfgangWaltenberger]

(if you allow me to chime in) Interesting discussion! For me, the words "expected limit" would have implied "the hypothesized limit computed before I make my observation, assuming I find no new physics", which is in contradiction to the word "a posteriori" which for me implies "after having made my observation". But I see in this context "expected limit" means something like "the limit, given my observation, under the assumption that what I measured is the expecation value of the SM background, thus under the assumption of the SM null hypothesis to hold". It's conceptually a weird thing, no? The limit of the "signal strength" parameter, given the data, under the assumption that that parameter is zero? I am not sure I know what I can do with this quantity, particularly in the context of searches, where I do not assume the SM hypothesis to hold. It's something like the sensitivity but trying to update my background estimates by assuming no new physics is in my data, which I cannot know. (Am I misunderstanding something?)

Anyways, even if I cannot fully wrap my head around its meaning, it's well-defined and we can use it to validate the results in our database, so thanks a lot for the clarification.