In [1]:
import ROOT
ROOT.RooMsgService.instance().setGlobalKillBelow(5)
%jsroot on

from matplotlib import pyplot
%matplotlib inline

Welcome to JupyROOT 6.07/07

Compute Limits using Hypothesis Test Inversion (CLs Limits).

In our last exercise we will compute an upper limit using the frequentist hypothesis test inversion method. We need to perform the hypothesis test for different parameter of interest points and compute the corresponding p-values. Since we are interesting in computing a limit, the test null hypothesis, that we want to disprove, is the in this case the S+B model, while the alternate hypothesis is the B only model. It is important to remember this, when we construct the hypothesis test calculator.

For doing the inversion RooStats provides an helper class, HypoTestInverter, which is constructed using one of the HypoTest calculator (Frequentist, Hybrid or Asymptotic). To compute the limit we need to do:

  • Create/read a workspace with S+B and B only model. If the B model does not exist we can create as in the previous exercise (setting the poi to zero).
  • Create an HypoTest calculator but using as null the S+B model and as alt the B only model
  • Configure the calculator if needed (set test statistics, number of toys, etc..)
  • Create the HypoTestInverter class and configure it (e.g. to set the number and range of points to scan, to use CLs for computing the limit, etc…)
  • Run the Inverter using GetInterval which returns an HypoTestInverterResult object.
  • Look at the result (e.g. limit, expected limits and bands)
  • Plot the result of the scan and the test statistic distribution obtained at every point.
In [2]:
infile = ROOT.TFile.Open("output/mH_BinModel.root")
w = infile.Get("w")
mc = w.obj("ModelConfig")
data = w.data("data")

sbModel = w.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()
bModel = sbModel.Clone()
bModel.SetName(sbModel.GetName()+"_with_poi_0")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
In [3]:
fc = ROOT.RooStats.FrequentistCalculator(data, bModel, sbModel)
fc.SetToys(1000,500)
ac = ROOT.RooStats.AsymptoticCalculator(data, bModel, sbModel)
ac.SetOneSided(True)
ac.SetPrintLevel(-1)

calc = ROOT.RooStats.HypoTestInverter(ac)
calc.SetConfidenceLevel(0.95)

AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1
 **********
 **    1 **SET PRINT           0
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 b            3.00000e+02  1.00000e+02    0.00000e+00  1.00000e+03
     2 s            0.00000e+00  5.00000e+01    0.00000e+00  5.00000e+02
 MINUIT WARNING IN PARAM DEF
 ============== STARTING VALUE IS AT LIMIT.
 MINUIT WARNING IN PARAMETR
 ============== VARIABLE2 IS AT ITS LOWER ALLOWED LIMIT.
 MINUIT WARNING IN PARAMETR
 ============== VARIABLE2 BROUGHT BACK INSIDE LIMITS.
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           0
 **********
 **********
 **    5 **SET STR           1
 **********
 **********
 **    6 **MIGRAD        1000           1
 **********
 MINUIT WARNING IN MIGrad    
 ============== VARIABLE2 IS AT ITS LOWER ALLOWED LIMIT.
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 FCN=8.96527 FROM MIGRAD    STATUS=CONVERGED      63 CALLS          64 TOTAL
                     EDM=8.14007e-06    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  b            3.01523e+02   6.01159e+01   9.00474e-05  -3.10875e-02
   2  s            1.24468e+02   6.28002e+01   2.04284e-04  -4.53302e-03
                               ERR DEF= 0.5
AsymptoticCalculator::EvaluateNLL -  value = 8.96527	fit time : Real time 0:00:00, CP time 0.100
MakeAsimov: Setting poi s to a constant value = 0
MakeAsimov: doing a conditional fit for finding best nuisance values 
 **********
 **    1 **SET PRINT           0
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 b            3.01523e+02  6.01159e+01    0.00000e+00  1.00000e+03
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           0
 **********
 **********
 **    5 **SET STR           1
 **********
 **********
 **    6 **MIGRAD         500           1
 **********
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 FCN=10.8691 FROM MIGRAD    STATUS=CONVERGED      18 CALLS          19 TOTAL
                     EDM=1.06662e-08    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  b            4.13325e+02   1.90015e+01   8.96260e-05  -2.67579e-03
                               ERR DEF= 0.5
fit time Real time 0:00:00, CP time 0.010
Generated Asimov data for observables RooArgSet:: = (nobs)
Generated Asimov data for global observables RooArgSet:: = (b0)
AsymptoticCalculator::EvaluateNLL  ........ using Minuit / Migrad with strategy  1 and tolerance 1
 **********
 **    7 **SET PRINT           0
 **********
 **********
 **    8 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 b            4.13325e+02  1.90015e+01    0.00000e+00  1.00000e+03
 **********
 **    9 **SET ERR         0.5
 **********
 **********
 **   10 **SET PRINT           0
 **********
 **********
 **   11 **SET STR           1
 **********
 **********
 **   12 **MIGRAD         500           1
 **********
 MIGRAD MINIMIZATION HAS CONVERGED.
 FCN=8.95017 FROM MIGRAD    STATUS=CONVERGED      15 CALLS          16 TOTAL
                     EDM=1.04963e-15    STRATEGY= 1  ERROR MATRIX UNCERTAINTY 100.0 per cent
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  b            4.13325e+02   3.59773e+00   6.04049e-11  -1.17107e-06
                               ERR DEF= 0.5
AsymptoticCalculator::EvaluateNLL -  value = 8.95017 for poi fixed at = 0	fit time : Real time 0:00:00, CP time 0.000

for CLS

In [4]:
calc.UseCLs(True)
calc.SetVerbose(False)
In [5]:
toymcs = calc.GetHypoTestCalculator().GetTestStatSampler()

Profile likelihood test statistics are used one sided for CLs since it's a bounded interval.

In [6]:
profll = ROOT.RooStats.ProfileLikelihoodTestStat(sbModel.GetPdf())
profll.SetOneSided(True)
toymcs.SetTestStatistic(profll)
In [7]:
if not sbModel.GetPdf().canBeExtended():
    toymcs.SetNEventsPerToy(1)
    print 'can not be extended'

can not be extended

define the number of points and some minmaxes

In [8]:
npoints = 10
poimin = poi.getMin()
poimax = poi.getMax()
In [9]:
calc.SetFixedScan(npoints,poimin,poimax)
r = calc.GetInterval()
upperLimit = r.UpperLimit()
In [10]:
plot = ROOT.RooStats.HypoTestInverterPlot("HTI_Result_Plot","HypoTest Scan Result",r)
c = ROOT.TCanvas("HypoTestInverter Scan")
c.SetLogy(False)
plot.Draw("CLb 2CL")
c.Draw()
In [ ]: