ROOT   6.21/01 Reference Guide
IntervalExamples.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// Example showing confidence intervals with four techniques.
5 ///
6 /// An example that shows confidence intervals with four techniques.
7 /// The model is a Normal Gaussian G(x|mu,sigma) with 100 samples of x.
8 /// The answer is known analytically, so this is a good example to validate
9 /// the RooStats tools.
10 ///
11 /// - expected interval is [-0.162917, 0.229075]
12 /// - plc interval is [-0.162917, 0.229075]
13 /// - fc interval is [-0.17 , 0.23] // stepsize is 0.01
14 /// - bc interval is [-0.162918, 0.229076]
15 /// - mcmc interval is [-0.166999, 0.230224]
16 ///
17 /// \macro_image
18 /// \macro_output
19 /// \macro_code
20 ///
21 /// \author Kyle Cranmer
22
23 #include "RooStats/ConfInterval.h"
32
33 #include "RooStats/ProofConfig.h"
34 #include "RooStats/ToyMCSampler.h"
35
36 #include "RooRandom.h"
37 #include "RooDataSet.h"
38 #include "RooRealVar.h"
39 #include "RooConstVar.h"
41 #include "RooDataHist.h"
42 #include "RooPoisson.h"
43 #include "RooPlot.h"
44
45 #include "TCanvas.h"
46 #include "TTree.h"
47 #include "TStyle.h"
48 #include "TMath.h"
49 #include "Math/DistFunc.h"
50 #include "TH1F.h"
51 #include "TMarker.h"
52 #include "TStopwatch.h"
53
54 #include <iostream>
55
57 using namespace RooFit;
58 using namespace RooStats;
59
60 void IntervalExamples()
61 {
62
63  // Time this macro
64  TStopwatch t;
65  t.Start();
66
67  // set RooFit random seed for reproducible results
69
70  // make a simple model via the workspace factory
71  RooWorkspace *wspace = new RooWorkspace();
72  wspace->factory("Gaussian::normal(x[-10,10],mu[-1,1],sigma[1])");
73  wspace->defineSet("poi", "mu");
74  wspace->defineSet("obs", "x");
75
76  // specify components of model for statistical tools
77  ModelConfig *modelConfig = new ModelConfig("Example G(x|mu,1)");
78  modelConfig->SetWorkspace(*wspace);
79  modelConfig->SetPdf(*wspace->pdf("normal"));
80  modelConfig->SetParametersOfInterest(*wspace->set("poi"));
81  modelConfig->SetObservables(*wspace->set("obs"));
82
83  // create a toy dataset
84  RooDataSet *data = wspace->pdf("normal")->generate(*wspace->set("obs"), 100);
85  data->Print();
86
87  // for convenience later on
88  RooRealVar *x = wspace->var("x");
89  RooRealVar *mu = wspace->var("mu");
90
91  // set confidence level
92  double confidenceLevel = 0.95;
93
94  // example use profile likelihood calculator
95  ProfileLikelihoodCalculator plc(*data, *modelConfig);
96  plc.SetConfidenceLevel(confidenceLevel);
97  LikelihoodInterval *plInt = plc.GetInterval();
98
99  // example use of Feldman-Cousins
100  FeldmanCousins fc(*data, *modelConfig);
101  fc.SetConfidenceLevel(confidenceLevel);
102  fc.SetNBins(100); // number of points to test per parameter
103  fc.UseAdaptiveSampling(true); // make it go faster
104
105  // Here, we consider only ensembles with 100 events
106  // The PDF could be extended and this could be removed
107  fc.FluctuateNumDataEntries(false);
108
109  // Proof
110  // ProofConfig pc(*wspace, 4, "workers=4", kFALSE); // proof-lite
111  // ProofConfig pc(w, 8, "localhost"); // proof cluster at "localhost"
112  // ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
113  // toymcsampler->SetProofConfig(&pc); // enable proof
114
115  PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();
116
117  // example use of BayesianCalculator
118  // now we also need to specify a prior in the ModelConfig
119  wspace->factory("Uniform::prior(mu)");
120  modelConfig->SetPriorPdf(*wspace->pdf("prior"));
121
122  // example usage of BayesianCalculator
123  BayesianCalculator bc(*data, *modelConfig);
124  bc.SetConfidenceLevel(confidenceLevel);
125  SimpleInterval *bcInt = bc.GetInterval();
126
127  // example use of MCMCInterval
128  MCMCCalculator mc(*data, *modelConfig);
129  mc.SetConfidenceLevel(confidenceLevel);
130  // special options
131  mc.SetNumBins(200); // bins used internally for representing posterior
132  mc.SetNumBurnInSteps(500); // first N steps to be ignored as burn-in
133  mc.SetNumIters(100000); // how long to run chain
134  mc.SetLeftSideTailFraction(0.5); // for central interval
135  MCMCInterval *mcInt = mc.GetInterval();
136
137  // for this example we know the expected intervals
138  double expectedLL =
139  data->mean(*x) + ROOT::Math::normal_quantile((1 - confidenceLevel) / 2, 1) / sqrt(data->numEntries());
140  double expectedUL =
141  data->mean(*x) + ROOT::Math::normal_quantile_c((1 - confidenceLevel) / 2, 1) / sqrt(data->numEntries());
142
143  // Use the intervals
144  std::cout << "expected interval is [" << expectedLL << ", " << expectedUL << "]" << endl;
145
146  cout << "plc interval is [" << plInt->LowerLimit(*mu) << ", " << plInt->UpperLimit(*mu) << "]" << endl;
147
148  std::cout << "fc interval is [" << interval->LowerLimit(*mu) << " , " << interval->UpperLimit(*mu) << "]" << endl;
149
150  cout << "bc interval is [" << bcInt->LowerLimit() << ", " << bcInt->UpperLimit() << "]" << endl;
151
152  cout << "mc interval is [" << mcInt->LowerLimit(*mu) << ", " << mcInt->UpperLimit(*mu) << "]" << endl;
153
154  mu->setVal(0);
155  cout << "is mu=0 in the interval? " << plInt->IsInInterval(RooArgSet(*mu)) << endl;
156
157  // make a reasonable style
164  gStyle->SetFillColor(0);
166  gStyle->SetStatColor(0);
167
168  // some plots
169  TCanvas *canvas = new TCanvas("canvas");
170  canvas->Divide(2, 2);
171
172  // plot the data
173  canvas->cd(1);
174  RooPlot *frame = x->frame();
175  data->plotOn(frame);
176  data->statOn(frame);
177  frame->Draw();
178
179  // plot the profile likelihood
180  canvas->cd(2);
181  LikelihoodIntervalPlot plot(plInt);
182  plot.Draw();
183
184  // plot the MCMC interval
185  canvas->cd(3);
186  MCMCIntervalPlot *mcPlot = new MCMCIntervalPlot(*mcInt);
187  mcPlot->SetLineColor(kGreen);
188  mcPlot->SetLineWidth(2);
189  mcPlot->Draw();
190
191  canvas->cd(4);
192  RooPlot *bcPlot = bc.GetPosteriorPlot();
193  bcPlot->Draw();
194
195  canvas->Update();
196
197  t.Stop();
198  t.Print();
199 }
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&)
Definition: RooAbsPdf.h:55
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Calls RooPlot* plotOn(RooPlot* frame, const RooLinkedList& cmdList) const ;.
Definition: RooAbsData.cxx:549
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void SetStatColor(Color_t color=19)
Definition: TStyle.h:371
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual void SetWorkspace(RooWorkspace &ws)
Definition: ModelConfig.h:66
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
Definition: TCanvas.cxx:696
This class provides simple and straightforward utilities to plot a LikelihoodInterval object...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
void Draw(const Option_t *options=NULL)
Definition: Rtypes.h:64
virtual void SetObservables(const RooArgSet &set)
Specify the observables.
Definition: ModelConfig.h:146
Double_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
void SetCanvasColor(Color_t color=19)
Definition: TStyle.h:325
double sqrt(double)
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
void SetFrameFillColor(Color_t color=1)
Definition: TStyle.h:353
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:81
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
virtual Double_t LowerLimit()
virtual Bool_t IsInInterval(const RooArgSet &) const
check if given point is in the interval
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:252
Double_t LowerLimit(RooRealVar &param)
return lower limit on a given parameter
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
Definition: TStyle.h:338
void SetCanvasBorderMode(Int_t mode=1)
Definition: TStyle.h:327
Double_t mean(const RooRealVar &var, const char *cutSpec=0, const char *cutRange=0) const
Definition: RooAbsData.h:193
Definition: TStyle.h:336
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual RooPlot * statOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Add a box with statistics information to the specified frame.
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
void SetLineColor(Color_t color)
void SetLineWidth(Int_t width)
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsData.h:175
PointSetInterval is a concrete implementation of the ConfInterval interface.
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:87
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
RooFactoryWSTool & factory()
Return instance to factory tool.
virtual Double_t UpperLimit()
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
void SetTitleFillColor(Color_t color=1)
Definition: TStyle.h:385
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
Definition: ModelConfig.h:100
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void Update()