ROOT   Reference Guide
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"
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
57using namespace RooFit;
58using namespace RooStats;
59
60void 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
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}
@ kGreen
Definition Rtypes.h:66
double sqrt(double)
R__EXTERN TStyle * gStyle
Definition TStyle.h:412
static struct mg_connection * fc(struct mg_context *ctx)
Definition civetweb.c:3728
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.
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition RooAbsData.h:193
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
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
Double_t mean(const RooRealVar &var, const char *cutSpec=0, const char *cutRange=0) const
Definition RooAbsData.h:211
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&,...
Definition RooAbsPdf.h:58
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
RooDataSet is a container class to hold unbinned data.
Definition RooDataSet.h:33
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:44
static RooPlot * frame(const RooAbsRealLValue &var, Double_t xmin, Double_t xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition RooPlot.cxx:249
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition RooPlot.cxx:691
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:53
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual void setVal(Double_t value)
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_t LowerLimit(const RooRealVar &param)
return the lower bound of the interval on a given parameter
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
virtual Bool_t IsInInterval(const RooArgSet &) const
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 SetLineWidth(Int_t width)
void Draw(const Option_t *options=NULL)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
virtual Double_t 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:30
virtual void SetObservables(const RooArgSet &set)
Specify the observables.
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the the workspace if not already there.
Definition ModelConfig.h:87
virtual void SetWorkspace(RooWorkspace &ws)
Definition ModelConfig.h:66
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition ModelConfig.h:81
PointSetInterval is a concrete implementation of the ConfInterval interface.
Double_t UpperLimit(RooRealVar &param)
return upper limit on a given parameter
Double_t 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_t UpperLimit()
virtual Double_t LowerLimit()
The RooWorkspace is a persistable container for RooFit projects.
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
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.
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...
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
The Canvas class.
Definition TCanvas.h:23
Definition TCanvas.cxx:708
void Update() override
Definition TCanvas.cxx:2504
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
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:608
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
Print the real and cpu time passed between the start and stop events.
Definition TStyle.h:340
void SetFrameFillColor(Color_t color=1)
Definition TStyle.h:355
void SetCanvasColor(Color_t color=19)
Definition TStyle.h:327
void SetCanvasBorderMode(Int_t mode=1)
Definition TStyle.h:329
void SetTitleFillColor(Color_t color=1)
Definition TStyle.h:387
void SetStatColor(Color_t color=19)
Definition TStyle.h:373