Limits from HypoTestResult::Significance() vs HypoTestInverter::GetInterval()

Roofit and RooStats
1

Dear Experts,

I’m trying to cross-check the result I get from HypoTestInverter with a 95% confidence interval. To do this, I am scanning over different significances for different values of the POI from HypoTestResult::Significance() until this goes over 1.64 (the standard deviation of a 5% single-sided Gaussian tail). I get different results from both of these, am I misunderstanding something?

(I’ve read through the code of HypoTestInverterResult and HypoTestInverter, which scan over pValue instead of significance, but I don’t see where the difference comes from)

I’ve attached a standalone script (demoSignificanceVsInverter.py (3.8 KB)). The key parts of this script are:

Strategy 1 (Using HypoTestResult::Significance())

# scan over values of mu until signifigance above 95% one-sided threshold
threshold = 1.64
for mu in np.arange(100,130,1):
    # set up model config
    h1_model = ModelConfig("h1",w)
    h1_model.SetPdf(w.pdf("model"))
    w.var("mu").setVal(mu)
    h1_model.SetObservables(RooArgSet(w.var("x")))
    h1_model.SetParametersOfInterest(w.set("poi"))
    h1_model.SetNuisanceParameters(w.set("nuis"))
    h1_model.SetSnapshot(w.set("poi"))
    # set up calculator
    ac = AsymptoticCalculator(data, h0_model, h1_model); # wrong order (for limits)
    ac.SetOneSided(True);
    result = ac.GetHypoTest()
    # if signifigance is higher, stop
    if result.Significance()>threshold:
        strategyResult1 = mu

Strategy 2 (using HypoTestInverter)

# set up calculator
ac = AsymptoticCalculator(data, h0_model, h1_model); # wrong order (for limits)
ac.SetOneSided(True);
calc = HypoTestInverter(ac);
calc.SetFixedScan(100,0,200, False)
calc.UseCLs(True)
calc.SetConfidenceLevel(0.95);
result = calc.GetInterval()
strategyResult2 = result.GetExpectedUpperLimit(0)

The issue is that strategyResult1 and strategyResult2 don’t agree:

==================================================
Limit on 'mu' from Strategy 1: 129
Limit on 'mu' from Strategy 2: 104.042149599
==================================================

Thanks for any insight about this difference.
Aaron

2

Hi @aaronsw,

for Strategy 2, you are using the CLs strategy. This re-normalises the p-values, and makes the limits more conservative. Check if that’s causing the differences.
Edit: I just realised that this would go in the other direction, unless you inverted mu.

3

Hi @StephanH,

Thanks! I changed to UseCLs(False) . Indeed it does make the limit stronger and farther away from “Strategy 1” :

==================================================
Limit on 'mu' from Strategy 1: 129
Limit on 'mu' from Strategy 2: 87.1695769378
==================================================

Could there be another explanation?
Thanks

4

I don’t know. Maybe @moneta can help?

5

Hi @moneta, do you have any ideas why I get different results with these two methods of estimating limits? I have a stand-alone illustration in the original post. Any help to understand what I’m doing wrong would be much appreciated.
Thanks,
Aaron

6

Hey Aaron,

I spoke to Lorenzo, and we noticed something:
From Strategy 1:

This retrieves the observed significance, the one that depends on the value of the test statistic in data.

From strategy 2:

This estimates an expected limit, which does not depend on data. Try UpperLimit() without CLs to have something that’s equivalent to case 1.

7

And here’s the fixed script:

################################################################################
# Strategy 1:
#   * Scan over different HypoTestResult::Significance() until reaching 1.64
################################################################################

# scan over values of mu until signifigance above 95% one-sided threshold
threshold = 1.6448
# set up model config
for mu in np.arange(100., 150. ,0.1):
    w.var("mu").setVal(mu)
    h1_model = ModelConfig("h1",w)
    h1_model.SetPdf(w.pdf("model"))
    h1_model.SetObservables(RooArgSet(w.var("x")))
    h1_model.SetParametersOfInterest(w.set("poi"))
    h1_model.SetNuisanceParameters(w.set("nuis"))
    h1_model.SetSnapshot(w.set("poi"))
    # set up calculator
    ac = AsymptoticCalculator(data, h0_model, h1_model); # wrong order (for limits)
    ac.SetOneSided(True);
    result = ac.GetHypoTest()
    # if signifigance is higher, stop
    print "Mu=", mu, "Signif=", result.Significance()
    if result.Significance()>threshold:
        strategyResult1 = mu
        break

################################################################################
# Strategy 2:
#   * Use hypotest inverter to get 95% CL limit
################################################################################

# set up h1 model
h1_model = ModelConfig("h1",w)
h1_model.SetPdf(w.pdf("model"))
w.var("mu").setVal(200)
h1_model.SetObservables(RooArgSet(w.var("x")))
h1_model.SetParametersOfInterest(w.set("poi"))
h1_model.SetNuisanceParameters(w.set("nuis"))
h1_model.SetSnapshot(w.set("poi"))

# set up calculator
ac = AsymptoticCalculator(data, h0_model, h1_model); # wrong order (for limits)
ac.SetOneSided(True);
calc = HypoTestInverter(ac);
calc.SetFixedScan(100,80,150, False)
calc.UseCLs(False)
calc.SetConfidenceLevel(0.95);
result = calc.GetInterval()
strategyResult2 = result.GetExpectedUpperLimit(0)


print "\n\n","="*50
print "Limit on 'mu' from Strategy 1:", strategyResult1
print "Limit on 'mu' from Strategy 2:", result.UpperLimit(), strategyResult2
print "="*50

Result:

==================================================
Limit on 'mu' from Strategy 1: 125.09999999999857
Limit on 'mu' from Strategy 2: 125.07723144 87.1625086103
==================================================
1 Like
8

Hi @StephanH (and Lorenzo), thank you so much for helping with this! This makes a lot of sense. I really appreciate it.
Best,
Aaron

9

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.