Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
FourBinInstructional.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roostats
3## \notebook
4## This example is a generalization of the on/off problem.
5##
6## This example is a generalization of the on/off problem.
7## It's a common setup for SUSY searches. Imagine that one has two
8## variables "x" and "y" (eg. missing ET and SumET), see figure.
9## The signal region has high values of both of these variables (top right).
10## One can see low values of "x" or "y" acting as side-bands. If we
11## just used "y" as a sideband, we would have the on/off problem.
12## - In the signal region we observe non events and expect s+b events.
13## - In the region with low values of "y" (bottom right)
14## we observe noff events and expect tau*b events.
15## Note the significance of tau. In the background only case:
16##
17## ~~~{.cpp}
18## tau ~ <expectation off> / <expectation on>
19## ~~~
20##
21## If tau is known, this model is sufficient, but often tau is not known exactly.
22## So one can use low values of "x" as an additional constraint for tau.
23## Note that this technique critically depends on the notion that the
24## joint distribution for "x" and "y" can be factorized.
25## Generally, these regions have many events, so it the ratio can be
26## measured very precisely there. So we extend the model to describe the
27## left two boxes... denoted with "bar".
28## - In the upper left we observe nonbar events and expect bbar events
29## - In the bottom left we observe noffbar events and expect tau bbar events
30## Note again we have:
31##
32## ~~~{.cpp}
33## tau ~ <expectation off bar> / <expectation on bar>
34## ~~~
35##
36## One can further expand the model to account for the systematic associated
37## to assuming the distribution of "x" and "y" factorizes (eg. that
38## tau is the same for off/on and offbar/onbar). This can be done in several
39## ways, but here we introduce an additional parameter rho, which so that
40## one set of models will use tau and the other tau*rho. The choice is arbitrary,
41## but it has consequences on the numerical stability of the algorithms.
42## The "bar" measurements typically have more events (& smaller relative errors).
43## If we choose
44##
45## ~~~{.cpp}
46## <expectation noffbar> = tau * rho * <expectation noonbar>
47## ~~~
48##
49## the product tau*rho will be known very precisely (~1/sqrt(bbar)) and the contour
50## in those parameters will be narrow and have a non-trivial tau~1/rho shape.
51## However, if we choose to put rho on the non/noff measurements (where the
52## product will have an error `~1/sqrt(b))`, the contours will be more amenable
53## to numerical techniques. Thus, here we choose to define:
54##
55## ~~~{.cpp}
56## tau := <expectation off bar> / (<expectation on bar>)
57## rho := <expectation off> / (<expectation on> * tau)
58##
59## ^ y
60## |
61## |---------------------------+
62## | | |
63## | nonbar | non |
64## | bbar | s+b |
65## | | |
66## |---------------+-----------|
67## | | |
68## | noffbar | noff |
69## | tau bbar | tau b rho |
70## | | |
71## +-----------------------------> x
72## ~~~
73##
74## Left in this way, the problem is under-constrained. However, one may
75## have some auxiliary measurement (usually based on Monte Carlo) to
76## constrain rho. Let us call this auxiliary measurement that gives
77## the nominal value of rho "rhonom". Thus, there is a 'constraint' term in
78## the full model: P(rhonom | rho). In this case, we consider a Gaussian
79## constraint with standard deviation sigma.
80##
81## In the example, the initial values of the parameters are:
82##
83## ~~~{.cpp}
84## - s = 40
85## - b = 100
86## - tau = 5
87## - bbar = 1000
88## - rho = 1
89## (sigma for rho = 20%)
90## ~~~
91##
92## and in the toy dataset:
93##
94## ~~~{.cpp}
95## - non = 139
96## - noff = 528
97## - nonbar = 993
98## - noffbar = 4906
99## - rhonom = 1.27824
100## ~~~
101##
102## Note, the covariance matrix of the parameters has large off-diagonal terms.
103## Clearly s,b are anti-correlated. Similarly, since noffbar >> nonbar, one would
104## expect bbar,tau to be anti-correlated.
105##
106## This can be seen below.
107##
108## ~~~{.cpp}
109## GLOBAL b bbar rho s tau
110## b 0.96820 1.000 0.191 -0.942 -0.762 -0.209
111## bbar 0.91191 0.191 1.000 0.000 -0.146 -0.912
112## rho 0.96348 -0.942 0.000 1.000 0.718 -0.000
113## s 0.76250 -0.762 -0.146 0.718 1.000 0.160
114## tau 0.92084 -0.209 -0.912 -0.000 0.160 1.000
115## ~~~
116##
117## Similarly, since tau*rho appears as a product, we expect rho,tau
118## to be anti-correlated. When the error on rho is significantly
119## larger than 1/sqrt(bbar), tau is essentially known and the
120## correlation is minimal (tau mainly cares about bbar, and rho about b,s).
121## In the alternate parametrization (bbar* tau * rho) the correlation coefficient
122## for rho,tau is large (and negative).
123##
124## The code below uses best-practices for RooFit & RooStats as of June 2010.
125##
126## It proceeds as follows:
127## - create a workspace to hold the model
128## - use workspace factory to quickly create the terms of the model
129## - use workspace factory to define total model (a prod pdf)
130## - create a RooStats ModelConfig to specify observables, parameters of interest
131## - add to the ModelConfig a prior on the parameters for Bayesian techniques
132## note, the pdf it is factorized for parameters of interest & nuisance params
133## - visualize the model
134## - write the workspace to a file
135## - use several of RooStats IntervalCalculators & compare results
136##
137## \macro_image
138## \macro_output
139## \macro_code
140##
141## \date July 2022
142## \authors Artem Busorgin, Kyle Cranmer and Tanja Rommerskirchen (C++ version)
143
144import ROOT
145
146doBayesian = False
147doFeldmanCousins = False
148doMCMC = False
149
150# let's time this challenging example
151t = ROOT.TStopwatch()
152t.Start()
153
154# set RooFit random seed for reproducible results
155ROOT.RooRandom.randomGenerator().SetSeed(4357)
156
157# make model
158wspace = ROOT.RooWorkspace("wspace")
159wspace.factory("Poisson::on(non[0,1000], sum::splusb(s[40,0,100],b[100,0,300]))")
160wspace.factory("Poisson::off(noff[0,5000], prod::taub(b,tau[5,3,7],rho[1,0,2]))")
161wspace.factory("Poisson::onbar(nonbar[0,10000], bbar[1000,500,2000])")
162wspace.factory("Poisson::offbar(noffbar[0,1000000], prod::lambdaoffbar(bbar, tau))")
163wspace.factory("Gaussian::mcCons(rhonom[1.,0,2], rho, sigma[.2])")
164wspace.factory("PROD::model(on,off,onbar,offbar,mcCons)")
165wspace.defineSet("obs", "non,noff,nonbar,noffbar,rhonom")
166
167wspace.factory("Uniform::prior_poi({s})")
168wspace.factory("Uniform::prior_nuis({b,bbar,tau, rho})")
169wspace.factory("PROD::prior(prior_poi,prior_nuis)")
170
171# ----------------------------------
172# Control some interesting variations
173# define parameers of interest
174# for 1-d plots
175wspace.defineSet("poi", "s")
176wspace.defineSet("nuis", "b,tau,rho,bbar")
177
178# for 2-d plots to inspect correlations:
179# wspace.defineSet("poi","s,rho")
180
181# test simpler cases where parameters are known.
182# wspace["tau"].setConstant()
183# wspace["rho"].setConstant()
184# wspace["b"].setConstant()
185# wspace["bbar"].setConstant()
186
187# inspect workspace
188# wspace.Print()
189
190# ----------------------------------------------------------
191# Generate toy data
192# generate toy data assuming current value of the parameters
193# import into workspace.
194# add Verbose() to see how it's being generated
195data = wspace["model"].generate(wspace.set("obs"), 1)
196# data.Print("v")
197wspace.Import(data)
198
199# ----------------------------------
200# Now the statistical tests
201# model config
202modelConfig = ROOT.RooStats.ModelConfig("FourBins")
203modelConfig.SetWorkspace(wspace)
204modelConfig.SetPdf(wspace["model"])
205modelConfig.SetPriorPdf(wspace["prior"])
206modelConfig.SetParametersOfInterest(wspace.set("poi"))
207modelConfig.SetNuisanceParameters(wspace.set("nuis"))
208wspace.Import(modelConfig)
209wspace.writeToFile("FourBin.root")
210
211# -------------------------------------------------
212# If you want to see the covariance matrix uncomment
213# wspace["model"].fitTo(data)
214
215# use ProfileLikelihood
216plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, modelConfig)
217plc.SetConfidenceLevel(0.95)
218plInt = plc.GetInterval()
219msglevel = ROOT.RooMsgService.instance().globalKillBelow()
220ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.FATAL)
221plInt.LowerLimit(wspace["s"]) # get ugly print out of the way. Fix.
222ROOT.RooMsgService.instance().setGlobalKillBelow(msglevel)
223
224# use FeldmaCousins (takes ~20 min)
225fc = ROOT.RooStats.FeldmanCousins(data, modelConfig)
226fc.SetConfidenceLevel(0.95)
227# number counting: dataset always has 1 entry with N events observed
228fc.FluctuateNumDataEntries(False)
229fc.UseAdaptiveSampling(True)
230fc.SetNBins(40)
231fcInt = ROOT.RooStats.PointSetInterval()
232if doFeldmanCousins: # takes 7 minutes
233 fcInt = fc.GetInterval()
234
235# use BayesianCalculator (only 1-d parameter of interest, slow for this problem)
236bc = ROOT.RooStats.BayesianCalculator(data, modelConfig)
237bc.SetConfidenceLevel(0.95)
238bInt = ROOT.RooStats.SimpleInterval()
239if doBayesian and len(wspace.set("poi")) == 1:
240 bInt = bc.GetInterval()
241else:
242 print("Bayesian Calc. only supports on parameter of interest")
243
244# use MCMCCalculator (takes about 1 min)
245# Want an efficient proposal function, so derive it from covariance
246# matrix of fit
247fit = wspace["model"].fitTo(data, Save=True)
248ph = ROOT.RooStats.ProposalHelper()
249ph.SetVariables(fit.floatParsFinal())
250ph.SetCovMatrix(fit.covarianceMatrix())
251ph.SetUpdateProposalParameters(ROOT.kTRUE) # auto-create mean vars and add mappings
252ph.SetCacheSize(100)
253pf = ph.GetProposalFunction()
254
255mc = ROOT.RooStats.MCMCCalculator(data, modelConfig)
256mc.SetConfidenceLevel(0.95)
257mc.SetProposalFunction(pf)
258mc.SetNumBurnInSteps(500) # first N steps to be ignored as burn-in
259mc.SetNumIters(50000)
260mc.SetLeftSideTailFraction(0.5) # make a central interval
261mcInt = ROOT.RooStats.MCMCInterval()
262if doMCMC:
263 mcInt = mc.GetInterval()
264
265# ----------------------------------
266# Make some plots
267c1 = ROOT.gROOT.Get("c1")
268if not c1:
269 c1 = ROOT.TCanvas("c1")
270
271if doBayesian and doMCMC:
272 c1.Divide(3)
273 c1.cd(1)
274elif doBayesian or doMCMC:
275 c1.Divide(2)
276 c1.cd(1)
277
278lrplot = ROOT.RooStats.LikelihoodIntervalPlot(plInt)
279lrplot.Draw()
280
281if doBayesian and len(wspace.set("poi")) == 1:
282 c1.cd(2)
283 # the plot takes a long time and print lots of error
284 # using a scan it is better
285 bc.SetScanOfPosterior(20)
286 bplot = bc.GetPosteriorPlot()
287 bplot.Draw()
288
289if doMCMC:
290 if doBayesian and len(wspace.set("poi")) == 1:
291 c1.cd(3)
292 else:
293 c1.cd(2)
294 mcPlot = ROOT.RooStats.MCMCIntervalPlot(mcInt)
295 mcPlot.Draw()
296
297# ----------------------------------
298# query intervals
299print(
300 "Profile Likelihood interval on s = [{}, {}]".format(plInt.LowerLimit(wspace["s"]), plInt.UpperLimit(wspace["s"]))
301)
302# Profile Likelihood interval on s = [12.1902, 88.6871]
303
304if doBayesian and len(wspace.set("poi")) == 1:
305 print("Bayesian interval on s = [{}, {}]".format(bInt.LowerLimit(), bInt.UpperLimit()))
306
307if doFeldmanCousins:
308 print(
309 "Feldman Cousins interval on s = [{}, {}]".format(fcInt.LowerLimit(wspace["s"]), fcInt.UpperLimit(wspace["s"]))
310 )
311 # Feldman Cousins interval on s = [18.75 +/- 2.45, 83.75 +/- 2.45]
312
313if doMCMC:
314 print("MCMC interval on s = [{}, {}]".format(mcInt.LowerLimit(wspace["s"]), mcInt.UpperLimit(wspace["s"])))
315 # MCMC interval on s = [15.7628, 84.7266]
316
317t.Stop()
318t.Print()
319
320c1.SaveAs("FourBinInstructional.png")
321
322# TODO: The calculators have to be destructed first. Otherwise, we can get
323# segmentation faults depending on the destruction order, which is random in
324# Python. Probably the issue is that some object has a non-owning pointer to
325# another object, which it uses in its destructor. This should be fixed either
326# in the design of RooStats in C++, or with phythonizations.
327del plc
328del bc
329del mc
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t UChar_t len
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t format