Logo ROOT  
Reference Guide
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook
4/// Example showing confidence intervals with four techniques.
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.
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]
17/// \macro_image
18/// \macro_output
19/// \macro_code
21/// \author Kyle Cranmer
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"
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"
54#include <iostream>
56// use this order for safety on library loading
57using namespace RooFit;
58using namespace RooStats;
60void IntervalExamples()
63 // Time this macro
64 TStopwatch t;
65 t.Start();
67 // set RooFit random seed for reproducible results
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");
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"));
83 // create a toy dataset
84 RooDataSet *data = wspace->pdf("normal")->generate(*wspace->set("obs"), 100);
85 data->Print();
87 // for convenience later on
88 RooRealVar *x = wspace->var("x");
89 RooRealVar *mu = wspace->var("mu");
91 // set confidence level
92 double confidenceLevel = 0.95;
94 // example use profile likelihood calculator
95 ProfileLikelihoodCalculator plc(*data, *modelConfig);
96 plc.SetConfidenceLevel(confidenceLevel);
97 LikelihoodInterval *plInt = plc.GetInterval();
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
105 // Here, we consider only ensembles with 100 events
106 // The PDF could be extended and this could be removed
107 fc.FluctuateNumDataEntries(false);
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
115 PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();
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"));
122 // example usage of BayesianCalculator
123 BayesianCalculator bc(*data, *modelConfig);
124 bc.SetConfidenceLevel(confidenceLevel);
125 SimpleInterval *bcInt = bc.GetInterval();
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();
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());
143 // Use the intervals
144 std::cout << "expected interval is [" << expectedLL << ", " << expectedUL << "]" << endl;
146 cout << "plc interval is [" << plInt->LowerLimit(*mu) << ", " << plInt->UpperLimit(*mu) << "]" << endl;
148 std::cout << "fc interval is [" << interval->LowerLimit(*mu) << " , " << interval->UpperLimit(*mu) << "]" << endl;
150 cout << "bc interval is [" << bcInt->LowerLimit() << ", " << bcInt->UpperLimit() << "]" << endl;
152 cout << "mc interval is [" << mcInt->LowerLimit(*mu) << ", " << mcInt->UpperLimit(*mu) << "]" << endl;
154 mu->setVal(0);
155 cout << "is mu=0 in the interval? " << plInt->IsInInterval(RooArgSet(*mu)) << endl;
157 // make a reasonable style
168 // some plots
169 TCanvas *canvas = new TCanvas("canvas");
170 canvas->Divide(2, 2);
172 // plot the data
173 canvas->cd(1);
174 RooPlot *frame = x->frame();
175 data->plotOn(frame);
176 data->statOn(frame);
177 frame->Draw();
179 // plot the profile likelihood
180 canvas->cd(2);
181 LikelihoodIntervalPlot plot(plInt);
182 plot.Draw();
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();
191 canvas->cd(4);
192 RooPlot *bcPlot = bc.GetPosteriorPlot();
193 bcPlot->Draw();
195 canvas->Update();
197 t.Stop();
198 t.Print();
@ kGreen
Definition: Rtypes.h:64
double sqrt(double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
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:175
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
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:546
Double_t mean(const RooRealVar &var, const char *cutSpec=0, const char *cutRange=0) const
Definition: RooAbsData.h:193
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:55
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
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:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:261
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.
Definition: MCMCInterval.h:30
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.
Definition: ModelConfig.h:146
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.
Definition: ModelConfig.h:100
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.
Definition: RooWorkspace.h:43
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:27
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2433
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:701
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)
Automatic pad generation by division.
Definition: TPad.cxx:1165
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
Stopwatch class.
Definition: TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
void SetPadBorderMode(Int_t mode=1)
Definition: TStyle.h:338
void SetFrameFillColor(Color_t color=1)
Definition: TStyle.h:353
void SetCanvasColor(Color_t color=19)
Definition: TStyle.h:325
void SetCanvasBorderMode(Int_t mode=1)
Definition: TStyle.h:327
void SetTitleFillColor(Color_t color=1)
Definition: TStyle.h:385
void SetStatColor(Color_t color=19)
Definition: TStyle.h:371
void SetPadColor(Color_t color=19)
Definition: TStyle.h:336
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
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19