Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
OneSidedFrequentistUpperLimitWithBands.py
Go to the documentation of this file.
1# \file
2# \ingroup tutorial_roostats
3# \notebook
4# OneSidedFrequentistUpperLimitWithBands
5#
6# This is a standard demo that can be used with any ROOT file
7# prepared in the standard way. You specify:
8# - name for input ROOT file
9# - name of workspace inside ROOT file that holds model and data
10# - name of ModelConfig that specifies details for calculator tools
11# - name of dataset
12#
13# With default parameters the macro will attempt to run the
14# standard hist2workspace example and read the ROOT file
15# that it produces.
16#
17# The first ~100 lines define a new test statistic, then the main macro starts.
18# You may want to control:
19# ~~~{.cpp}
20# double confidenceLevel=0.95;
21# int nPointsToScan = 12;
22# int nToyMC = 150;
23# ~~~
24# This uses a modified version of the profile likelihood ratio as
25# a test statistic for upper limits (eg. test stat = 0 if muhat>mu).
26#
27# Based on the observed data, one defines a set of parameter points
28# to be tested based on the value of the parameter of interest
29# and the conditional MLE (eg. profiled) values of the nuisance parameters.
30#
31# At each parameter point, pseudo-experiments are generated using this
32# fixed reference model and then the test statistic is evaluated.
33# Note, the nuisance parameters are floating in the fits. For each point,
34# the threshold that defines the 95% acceptance region is found. This
35# forms a "Confidence Belt".
36#
37# After constructing the confidence belt, one can find the confidence
38# interval for any particular dataset by finding the intersection
39# of the observed test statistic and the confidence belt. First
40# this is done on the observed data to get an observed 1-sided upper limt.
41#
42# Finally, there expected limit and bands (from background-only) are
43# formed by generating background-only data and finding the upper limit.
44# This is done by hand for now, will later be part of the RooStats tools.
45#
46# On a technical note, this technique is NOT the Feldman-Cousins technique,
47# because that is a 2-sided interval BY DEFINITION. However, like the
48# Feldman-Cousins technique this is a Neyman-Construction. For technical
49# reasons the easiest way to implement this right now is to use the
50# FeldmanCousins tool and then change the test statistic that it is using.
51#
52# Building the confidence belt can be computationally expensive. Once it is built,
53# one could save it to a file and use it in a separate step.
54#
55# Note, if you have a boundary on the parameter of interest (eg. cross-section)
56# the threshold on the one-sided test statistic starts off very small because we
57# are only including downward fluctuations. You can see the threshold in these printouts:
58# ~~~{.cpp}
59# [#0] PROGRESS:Generation -- generated toys: 500 / 999
60# NeymanConstruction: Prog: 12/50 total MC = 39 this test stat = 0
61# SigXsecOverSM=0.69 alpha_syst1=0.136515 alpha_syst3=0.425415 beta_syst2=1.08496 [-1e+30, 0.011215] in interval = 1
62# ~~~
63# this tells you the values of the parameters being used to generate the pseudo-experiments
64# and the threshold in this case is 0.011215. One would expect for 95% that the threshold
65# would be ~1.35 once the cross-section is far enough away from 0 that it is essentially
66# unaffected by the boundary. As one reaches the last points in the scan, the
67# threshold starts to get artificially high. This is because the range of the parameter in
68# the fit is the same as the range in the scan. In the future, these should be independently
69# controlled, but they are not now. As a result the ~50% of pseudo-experiments that have an
70# upward fluctuation end up with muhat = muMax. Because of this, the upper range of the
71# parameter should be well above the expected upper limit... but not too high or one will
72# need a very large value of nPointsToScan to resolve the relevant region. This can be
73# improved, but this is the first version of this script.
74#
75# Important note: when the model includes external constraint terms, like a Gaussian
76# constraint to a nuisance parameter centered around some nominal value there is
77# a subtlety. The asymptotic results are all based on the assumption that all the
78# measurements fluctuate... including the nominal values from auxiliary measurements.
79# If these do not fluctuate, this corresponds to an "conditional ensemble". The
80# result is that the distribution of the test statistic can become very non-chi^2.
81# This results in thresholds that become very large. This can be seen in the following
82# thought experiment. Say the model is
83# \f$ Pois(N | s + b)G(b0|b,sigma) \f$
84# where \f$ G(b0|b,sigma) \f$ is the external constraint and b0 is 100. If N is also 100
85# then the profiled value of b given s is going to be some trade off between 100-s and b0.
86# If sigma is \f$ \sqrt(N) \f$, then the profiled value of b is probably 100 - s/2 So for
87# s=60 we are going to have a profiled value of b~70. Now when we generate pseudo-experiments
88# for s=60, b=70 we will have N~130 and the average shat will be 30, not 60. In practice,
89# this is only an issue for values of s that are very excluded. For values of s near the 95%
90# limit this should not be a big effect. This can be avoided if the nominal values of the constraints also fluctuate,
91# but that requires that those parameters are RooRealVars in the model.
92# This version does not deal with this issue, but it will be addressed in a future version.
93#
94# \macro_image
95# \macro_output
96# \macro_code
97#
98# \authors Kyle Cranmer (C++ version), Haichen Wang, Daniel Whiteson, and P. P. (Python translation)
99
100import ROOT
101
102# -------------------------------------------------------
103# The actual macro
104
105
107 infile="", workspaceName="combined", modelConfigName="ModelConfig", dataName="obsData"
108):
109
110 confidenceLevel = 0.95
111 nPointsToScan = 12
112 nToyMC = 150
113
114 # -------------------------------------------------------
115 # First part is just to access a user-defined file
116 # or create the standard example file if it doesn't exist
117 filename = ""
118 if infile == "":
119 filename = "results/example_combined_GaussExample_model.root"
120 fileExist = not ROOT.gSystem.AccessPathName(filename) # note opposite return code
121 # if file does not exists generate with histfactory
122 if not fileExist:
123 # Normally this would be run on the command line
124 print(f"will run standard hist2workspace example")
125 ROOT.gROOT.ProcessLine(".not prepareHistFactory .")
126 ROOT.gROOT.ProcessLine(".not hist2workspace config/example.xml")
127 print(f"\n\n---------------------")
128 print(f"Done creating example input")
129 print(f"---------------------\n\n")
130
131 else:
132 filename = infile
133
134 # Try to open the file
135 file = ROOT.TFile.Open(filename)
136
137 # if input file was specified but not found, quit
138 if not file:
139 print(f"StandardRooStatsDemoMacro: Input file {filename} is not found")
140 return
141
142 # -------------------------------------------------------
143 # Now get the data and workspace
144
145 # get the workspace out of the file
146 w = file.Get(workspaceName)
147 global gw
148 gw = w
149 global gfile
150 gfile = file
151
152 if not w:
153 print(f"workspace not found")
154 return
155
156 # get the modelConfig out of the file
157 mc = w.obj(modelConfigName)
158
159 # get the modelConfig out of the file
160 data = w.data(dataName)
161
162 # make sure ingredients are found
163 if not data or not mc:
164 w.Print()
165 print(f"data or ModelConfig was not found")
166 return
167
168 # -------------------------------------------------------
169 # Now get the POI for convenience
170 # you may want to adjust the range of your POI
171
172 firstPOI = mc.GetParametersOfInterest().first()
173 # firstPOI->setMin(0);
174 # firstPOI->setMax(10);
175
176 # --------------------------------------------
177 # Create and use the FeldmanCousins tool
178 # to find and plot the 95% confidence interval
179 # on the parameter of interest as specified
180 # in the model config
181 # REMEMBER, we will change the test statistic
182 # so this is NOT a Feldman-Cousins interval
183 fc = ROOT.RooStats.FeldmanCousins(data, mc)
184 fc.SetConfidenceLevel(confidenceLevel)
185 fc.AdditionalNToysFactor(
186 0.5
187 ) # degrade/improve sampling that defines confidence belt: in this case makes the example faster
188 # fc.UseAdaptiveSampling(True); # speed it up a bit, don't use for expected limits
189 fc.SetNBins(nPointsToScan) # set how many points per parameter of interest to scan
190 fc.CreateConfBelt(True) # save the information in the belt for plotting
191
192 # -------------------------------------------------------
193 # Feldman-Cousins is a unified limit by definition
194 # but the tool takes care of a few things for us like which values
195 # of the nuisance parameters should be used to generate toys.
196 # so let's just change the test statistic and realize this is
197 # no longer "Feldman-Cousins" but is a fully frequentist Neyman-Construction.
198 # ProfileLikelihoodTestStatModified onesided(*mc->GetPdf());
199 # fc.GetTestStatSampler()->SetTestStatistic(&onesided);
200 # ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(True);
201 toymcsampler = fc.GetTestStatSampler()
202 testStat = toymcsampler.GetTestStatistic()
203 testStat.SetOneSided(True)
204
205 # Since this tool needs to throw toy MC the PDF needs to be
206 # extended or the tool needs to know how many entries in a dataset
207 # per pseudo experiment.
208 # In the 'number counting form' where the entries in the dataset
209 # are counts, and not values of discriminating variables, the
210 # datasets typically only have one entry and the PDF is not
211 # extended.
212 if not mc.GetPdf().canBeExtended():
213 if data.numEntries() == 1:
214 fc.FluctuateNumDataEntries(False)
215 else:
216 print(f"Not sure what to do about this model")
217
218 if mc.GetGlobalObservables():
219 print(f"will use global observables for unconditional ensemble")
220 mc.GetGlobalObservables().Print()
221 toymcsampler.SetGlobalObservables(mc.GetGlobalObservables())
222
223 # Now get the interval
224 interval = fc.GetInterval()
225 belt = fc.GetConfidenceBelt()
226
227 # print out the interval on the first Parameter of Interest
228 print(
229 f"\n95% interval on {firstPOI.GetName()} is : [{interval.LowerLimit(firstPOI)}, {interval.UpperLimit(firstPOI)} ]"
230 )
231
232 # get observed UL and value of test statistic evaluated there
233 tmpPOI = ROOT.RooArgSet(firstPOI)
234 observedUL = interval.UpperLimit(firstPOI)
235 firstPOI.setVal(observedUL)
236 obsTSatObsUL = fc.GetTestStatSampler().EvaluateTestStatistic(data, tmpPOI)
237
238 # Ask the calculator which points were scanned
239 parameterScan = fc.GetPointsToScan()
240 tmpPoint = ROOT.RooArgSet()
241
242 # make a histogram of parameter vs. threshold
243 histOfThresholds = ROOT.TH1F(
244 "histOfThresholds", "", parameterScan.numEntries(), firstPOI.getMin(), firstPOI.getMax()
245 )
246 histOfThresholds.GetXaxis().SetTitle(firstPOI.GetName())
247 histOfThresholds.GetYaxis().SetTitle("Threshold")
248
249 # loop through the points that were tested and ask confidence belt
250 # what the upper/lower thresholds were.
251 # For FeldmanCousins, the lower cut off is always 0
252 for i in range(parameterScan.numEntries()):
253 tmpPoint = parameterScan.get(i).clone("temp")
254 # cout <<"get threshold"<<endl;
255 arMax = belt.GetAcceptanceRegionMax(tmpPoint)
256 poiVal = tmpPoint.getRealValue(firstPOI.GetName())
257 histOfThresholds.Fill(poiVal, arMax)
258
259 c1 = ROOT.TCanvas()
260 c1.Divide(2)
261 c1.cd(1)
262 histOfThresholds.SetMinimum(0)
263 histOfThresholds.Draw()
264 c1.Update()
265 c1.Draw()
266 c1.cd(2)
267
268 # -------------------------------------------------------
269 # Now we generate the expected bands and power-constraint
270
271 # First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
272 nll = mc.GetPdf().createNLL(data)
273 profile = nll.createProfile(mc.GetParametersOfInterest())
274 firstPOI.setVal(0.0)
275 profile.getVal() # this will do fit and set nuisance parameters to profiled values
276 poiAndNuisance = ROOT.RooArgSet()
277 if mc.GetNuisanceParameters():
278 poiAndNuisance.add(mc.GetNuisanceParameters())
279 poiAndNuisance.add(mc.GetParametersOfInterest())
280 w.saveSnapshot("paramsToGenerateData", poiAndNuisance)
281 paramsToGenerateData = poiAndNuisance.snapshot()
282 print(f"\nWill use these parameter points to generate pseudo data for bkg only")
283 paramsToGenerateData.Print("v")
284
285 unconditionalObs = ROOT.RooArgSet()
286 unconditionalObs.add(mc.GetObservables())
287 unconditionalObs.add(mc.GetGlobalObservables()) # comment this out for the original conditional ensemble
288
289 CLb = 0
290 CLbinclusive = 0
291
292 # Now we generate background only and find distribution of upper limits
293 histOfUL = ROOT.TH1F("histOfUL", "", 100, 0, firstPOI.getMax())
294 histOfUL.GetXaxis().SetTitle("Upper Limit (background only)")
295 histOfUL.GetYaxis().SetTitle("Entries")
296 for imc in range(nToyMC):
297
298 # set parameters back to values for generating pseudo data
299 # cout << "\n get current nuis, set vals, print again" << endl;
300 w.loadSnapshot("paramsToGenerateData")
301 # poiAndNuisance->Print("v");
302
303 toyData = ROOT.RooDataSet()
304 # debugging
305 global gmc
306 gmc = mc
307 # return
308 # now generate a toy dataset
309 if not mc.GetPdf().canBeExtended():
310 if data.numEntries() == 1:
311 toyData = mc.GetPdf().generate(mc.GetObservables(), 1)
312 else:
313 print(f"Not sure what to do about this model")
314 else:
315 # print("generating extended dataset")
316 toyData = mc.GetPdf().generate(mc.GetObservables(), Extended=True)
317
318 # generate global observables
319 # need to be careful for simpdf
320 # RooDataSet* globalData = mc->GetPdf()->generate(*mc->GetGlobalObservables(),1);
321
322 simPdf = mc.GetPdf()
323 if not simPdf:
324 one = mc.GetPdf().generate(mc.GetGlobalObservables(), 1)
325 values = one.get()
326 allVars = mc.GetPdf().getVariables()
327 allVars.assign(values)
328 # del values
329 # del one
330 else:
331
332 # try fix for sim pdf
333 for tt in simPdf.indexCat():
334 catName = tt.first
335 # global gcatName
336 # gcatName = catName
337 # return
338 # Get pdf associated with state from simpdf
339 pdftmp = simPdf.getPdf(str(catName))
340
341 # Generate only global variables defined by the pdf associated with this state
342 globtmp = pdftmp.getObservables(mc.GetGlobalObservables())
343 tmp = pdftmp.generate(globtmp, 1)
344
345 # Transfer values to output placeholder
346 globtmp.assign(tmp.get(0))
347
348 # globalData->Print("v");
349 # unconditionalObs = *globalData->get();
350 # mc->GetGlobalObservables()->Print("v");
351 # delete globalData;
352 # cout << "toy data = " << endl;
353 # toyData->get()->Print("v");
354
355 # get test stat at observed UL in observed data
356 firstPOI.setVal(observedUL)
357 toyTSatObsUL = fc.GetTestStatSampler().EvaluateTestStatistic(toyData, tmpPOI)
358 # toyData->get()->Print("v");
359 # cout <<"obsTSatObsUL " <<obsTSatObsUL << "toyTS " << toyTSatObsUL << endl;
360 if obsTSatObsUL < toyTSatObsUL: # not sure about <= part yet
361 CLb += (1.0) / nToyMC
362 if obsTSatObsUL <= toyTSatObsUL: # not sure about <= part yet
363 CLbinclusive += (1.0) / nToyMC
364
365 # loop over points in belt to find upper limit for this toy data
366 thisUL = ROOT.Double_t(0)
367 for i in range(parameterScan.numEntries()):
368 tmpPoint = parameterScan.get(i).clone("temp")
369 arMax = belt.GetAcceptanceRegionMax(tmpPoint)
370 firstPOI.setVal(tmpPoint.getRealValue(firstPOI.GetName()))
371 # double thisTS = profile->getVal();
372 thisTS = fc.GetTestStatSampler().EvaluateTestStatistic(toyData, tmpPOI)
373
374 # cout << "poi = " << firstPOI->getVal()
375 # << " max is " << arMax << " this profile = " << thisTS << endl;
376 # cout << "thisTS = " << thisTS<<endl;
377 if thisTS <= arMax:
378 thisUL = firstPOI.getVal()
379 else:
380 break
381
382 """
383 #
384 # loop over points in belt to find upper limit for this toy data
385 thisUL = 0
386 for i in range(histOfThresholds.GetNbinsX() ++i)
387 tmpPoint = (RooArgSet) parameterScan.get(i).clone("temp")
388 print("---------------- ", i)
389 tmpPoint.Print("v")
390 print(f"from hist ", histOfThresholds.GetBinCenter(i+1) )
391 arMax = histOfThresholds.GetBinContent(i+1)
392 # cout << " threshold from Hist = aMax " << arMax<<endl;
393 # double arMax2 = belt->GetAcceptanceRegionMax(*tmpPoint);
394 # cout << "from scan arMax2 = "<< arMax2 << endl; # not the same due to TH1F not TH1D
395 # cout << "scan - hist" << arMax2-arMax << endl;
396 firstPOI.setVal( histOfThresholds.GetBinCenter(i+1))
397 # double thisTS = profile->getVal();
398 thisTS = fc.GetTestStatSampler().EvaluateTestStatistic(toyData,tmpPOI)
399
400 # cout << "poi = " << firstPOI->getVal()
401 # = ROOT.Double_t() << " max is " << arMax << " this profile = " << thisTS << endl;
402 # cout << "thisTS = " << thisTS<<endl;
403
404 # NOTE: need to add a small epsilon term for single precision vs. double precision
405# if(thisTS<=arMax + 1e-7){
406# thisUL = firstPOI->getVal();
407# } else{
408# break;
409# }
410# }
411# */
412#
413 """
414
415 histOfUL.Fill(thisUL)
416
417 # for few events, data is often the same, and UL is often the same
418 # cout << "thisUL = " << thisUL<<endl;
419
420 # delete toyData
421 c1.cd(2)
422 histOfUL.Draw()
423 c1.Update()
424 c1.Draw()
425 c1.SaveAs("OneSidedFrequentistUpperLimitWithBands.png")
426
427 # if you want to see a plot of the sampling distribution for a particular scan point:
428 #
429 """
430 SamplingDistPlot sampPlot
431 indexInScan = 0
432 tmpPoint = (RooArgSet) parameterScan.get(indexInScan).clone("temp")
433 firstPOI.setVal( tmpPoint.getRealValue(firstPOI.GetName()) )
434 toymcsampler.SetParametersForTestStat(tmpPOI)
435 samp = toymcsampler.GetSamplingDistribution(tmpPoint)
436 sampPlot.AddSamplingDistribution(samp)
437 sampPlot.Draw()
438 """
439
440 # Now find bands and power constraint
441 bins = histOfUL.GetIntegral()
442 cumulative = histOfUL.Clone("cumulative")
443 cumulative.SetContent(bins)
444 band2sigDown = band1sigDown = bandMedian = band1sigUp = band2sigUp = ROOT.Double_t()
445 for i in range(cumulative.GetNbinsX()):
446 if bins[i] < ROOT.RooStats.SignificanceToPValue(2):
447 band2sigDown = cumulative.GetBinCenter(i)
448 if bins[i] < ROOT.RooStats.SignificanceToPValue(1):
449 band1sigDown = cumulative.GetBinCenter(i)
450 if bins[i] < 0.5:
451 bandMedian = cumulative.GetBinCenter(i)
452 if bins[i] < ROOT.RooStats.SignificanceToPValue(-1):
453 band1sigUp = cumulative.GetBinCenter(i)
454 if bins[i] < ROOT.RooStats.SignificanceToPValue(-2):
455 band2sigUp = cumulative.GetBinCenter(i)
456
457 print(f"-2 sigma band ", band2sigDown)
458 print(f"-1 sigma band {band1sigDown} [Power Constraint)]")
459 print(f"median of band ", bandMedian)
460 print(f"+1 sigma band ", band1sigUp)
461 print(f"+2 sigma band ", band2sigUp)
462
463 # print out the interval on the first Parameter of Interest
464 print(f"\nObserved 95% upper-limit ", interval.UpperLimit(firstPOI))
465 print(f"CLb strict [P(toy>obs|0)] for observed 95% upper-limit ", CLb)
466 print("inclusive [P(toy>=obs|0)] for observed 95% upper-limit ", CLbinclusive)
467
468
470 infile="", workspaceName="combined", modelConfigName="ModelConfig", dataName="obsData"
471)
void Print(GNN_Data &d, std::string txt="")