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

Roofit and RooStats
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