Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
32
35
36#include "RooRandom.h"
37#include "RooDataSet.h"
38#include "RooRealVar.h"
39#include "RooConstVar.h"
40#include "RooAddition.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
56// use this order for safety on library loading
57using namespace RooFit;
58using namespace RooStats;
59
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 std::unique_ptr<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
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);
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}
@ kGreen
Definition Rtypes.h:66
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
RooFit::OwningPtr< RooDataSet > generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={})
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition RooAbsPdf.h:57
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
static RooPlot * frame(const RooAbsRealLValue &var, double xmin, double xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:237
void Draw(Option_t *options=nullptr) override
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:649
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:48
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
The FeldmanCousins class (like the Feldman-Cousins technique) is essentially a specific configuration...
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
double UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
double LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
bool IsInInterval(const RooArgSet &) const override
check if given point is in the interval
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void SetLineColor(Color_t color)
void Draw(const Option_t *options=nullptr) override
void SetLineWidth(Int_t width)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual double UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
virtual double LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
virtual void SetObservables(const RooArgSet &set)
Specify the observables.
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the workspace if not already there.
Definition ModelConfig.h:95
virtual void SetWorkspace(RooWorkspace &ws)
Definition ModelConfig.h:71
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the workspace if not already there.
Definition ModelConfig.h:88
PointSetInterval is a concrete implementation of the ConfInterval interface.
double UpperLimit(RooRealVar &param)
return upper limit on a given parameter
double LowerLimit(RooRealVar &param)
return lower limit on a given parameter
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual double UpperLimit()
return the interval upper limit
virtual double LowerLimit()
return the interval lower limit
Persistable container for RooFit projects.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
const RooArgSet * set(RooStringView name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
RooFactoryWSTool & factory()
Return instance to factory tool.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
bool defineSet(const char *name, const RooArgSet &aset, bool importMissing=false)
Define a named RooArgSet with given constituents.
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
The Canvas class.
Definition TCanvas.h:23
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:716
void Update() override
Update canvas pad buffers.
Definition TCanvas.cxx:2483
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1153
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:615
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Stop()
Stop the stopwatch.
void Print(Option_t *option="") const override
Print the real and cpu time passed between the start and stop events.
void SetPadBorderMode(Int_t mode=1)
Definition TStyle.h:354
void SetFrameFillColor(Color_t color=1)
Definition TStyle.h:369
void SetCanvasColor(Color_t color=19)
Definition TStyle.h:341
void SetCanvasBorderMode(Int_t mode=1)
Definition TStyle.h:343
void SetTitleFillColor(Color_t color=1)
Definition TStyle.h:401
void SetStatColor(Color_t color=19)
Definition TStyle.h:387
void SetPadColor(Color_t color=19)
Definition TStyle.h:352
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
Double_t x[n]
Definition legend1.C:17
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
Namespace for the RooStats classes.
Definition Asimov.h:19