Logo ROOT  
Reference Guide
StandardBayesianMCMCDemo.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_roostats
3/// \notebook -js
4/// Standard demo of the Bayesian MCMC calculator
5///
6/// This is a standard demo that can be used with any ROOT file
7/// prepared in the standard way. You specify:
8/// - name for input ROOT file
9/// - name of workspace inside ROOT file that holds model and data
10/// - name of ModelConfig that specifies details for calculator tools
11/// - name of dataset
12///
13/// With default parameters the macro will attempt to run the
14/// standard hist2workspace example and read the ROOT file
15/// that it produces.
16///
17/// The actual heart of the demo is only about 10 lines long.
18///
19/// The MCMCCalculator is a Bayesian tool that uses
20/// the Metropolis-Hastings algorithm to efficiently integrate
21/// in many dimensions. It is not as accurate as the BayesianCalculator
22/// for simple problems, but it scales to much more complicated cases.
23///
24/// \macro_image
25/// \macro_output
26/// \macro_code
27///
28/// \author Kyle Cranmer
29
30#include "TFile.h"
31#include "TROOT.h"
32#include "TCanvas.h"
33#include "TMath.h"
34#include "TSystem.h"
35#include "RooWorkspace.h"
36#include "RooAbsData.h"
37
45#include "RooFitResult.h"
46
47using namespace RooFit;
48using namespace RooStats;
49
50struct BayesianMCMCOptions {
51
52 double confLevel = 0.95;
53 int intervalType = 2; // type of interval (0 is shortest, 1 central, 2 upper limit)
54 double maxPOI = -999; // force different values of POI for doing the scan (default is given value)
55 double minPOI = -999;
56 int numIters = 100000; // number of iterations
57 int numBurnInSteps = 100; // number of burn in steps to be ignored
58};
59
60BayesianMCMCOptions optMCMC;
61
62void StandardBayesianMCMCDemo(const char *infile = "", const char *workspaceName = "combined",
63 const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
64{
65
66 // -------------------------------------------------------
67 // First part is just to access a user-defined file
68 // or create the standard example file if it doesn't exist
69
70 const char *filename = "";
71 if (!strcmp(infile, "")) {
72 filename = "results/example_combined_GaussExample_model.root";
73 bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
74 // if file does not exists generate with histfactory
75 if (!fileExist) {
76#ifdef _WIN32
77 cout << "HistFactory file cannot be generated on Windows - exit" << endl;
78 return;
79#endif
80 // Normally this would be run on the command line
81 cout << "will run standard hist2workspace example" << endl;
82 gROOT->ProcessLine(".! prepareHistFactory .");
83 gROOT->ProcessLine(".! hist2workspace config/example.xml");
84 cout << "\n\n---------------------" << endl;
85 cout << "Done creating example input" << endl;
86 cout << "---------------------\n\n" << endl;
87 }
88
89 } else
90 filename = infile;
91
92 // Try to open the file
93 TFile *file = TFile::Open(filename);
94
95 // if input file was specified byt not found, quit
96 if (!file) {
97 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
98 return;
99 }
100
101 // -------------------------------------------------------
102 // Tutorial starts here
103 // -------------------------------------------------------
104
105 // get the workspace out of the file
106 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
107 if (!w) {
108 cout << "workspace not found" << endl;
109 return;
110 }
111
112 // get the modelConfig out of the file
113 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
114
115 // get the modelConfig out of the file
116 RooAbsData *data = w->data(dataName);
117
118 // make sure ingredients are found
119 if (!data || !mc) {
120 w->Print();
121 cout << "data or ModelConfig was not found" << endl;
122 return;
123 }
124
125 // Want an efficient proposal function
126 // default is uniform.
127
128 /*
129 // this one is based on the covariance matrix of fit
130 RooFitResult* fit = mc->GetPdf()->fitTo(*data,Save());
131 ProposalHelper ph;
132 ph.SetVariables((RooArgSet&)fit->floatParsFinal());
133 ph.SetCovMatrix(fit->covarianceMatrix());
134 ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
135 ph.SetCacheSize(100);
136 ProposalFunction* pf = ph.GetProposalFunction();
137 */
138
139 // this proposal function seems fairly robust
140 SequentialProposal sp(0.1);
141 // -------------------------------------------------------
142 // create and use the MCMCCalculator
143 // to find and plot the 95% credible interval
144 // on the parameter of interest as specified
145 // in the model config
146 MCMCCalculator mcmc(*data, *mc);
147 mcmc.SetConfidenceLevel(optMCMC.confLevel); // 95% interval
148 // mcmc.SetProposalFunction(*pf);
149 mcmc.SetProposalFunction(sp);
150 mcmc.SetNumIters(optMCMC.numIters); // Metropolis-Hastings algorithm iterations
151 mcmc.SetNumBurnInSteps(optMCMC.numBurnInSteps); // first N steps to be ignored as burn-in
152
153 // default is the shortest interval.
154 if (optMCMC.intervalType == 0)
155 mcmc.SetIntervalType(MCMCInterval::kShortest); // for shortest interval (not really needed)
156 if (optMCMC.intervalType == 1)
157 mcmc.SetLeftSideTailFraction(0.5); // for central interval
158 if (optMCMC.intervalType == 2)
159 mcmc.SetLeftSideTailFraction(0.); // for upper limit
160
161 RooRealVar *firstPOI = (RooRealVar *)mc->GetParametersOfInterest()->first();
162 if (optMCMC.minPOI != -999)
163 firstPOI->setMin(optMCMC.minPOI);
164 if (optMCMC.maxPOI != -999)
165 firstPOI->setMax(optMCMC.maxPOI);
166
167 MCMCInterval *interval = mcmc.GetInterval();
168
169 // make a plot
170 // TCanvas* c1 =
171 auto c1 = new TCanvas("IntervalPlot");
172 MCMCIntervalPlot plot(*interval);
173 plot.Draw();
174
175 TCanvas *c2 = new TCanvas("extraPlots");
176 const RooArgSet *list = mc->GetNuisanceParameters();
177 if (list->getSize() > 1) {
178 double n = list->getSize();
179 int ny = TMath::CeilNint(sqrt(n));
180 int nx = TMath::CeilNint(double(n) / ny);
181 c2->Divide(nx, ny);
182 }
183
184 // draw a scatter plot of chain results for poi vs each nuisance parameters
186 RooRealVar *nuis = NULL;
187 int iPad = 1; // iPad, that's funny
188 while ((nuis = (RooRealVar *)it->Next())) {
189 c2->cd(iPad++);
190 plot.DrawChainScatter(*firstPOI, *nuis);
191 }
192
193 // print out the interval on the first Parameter of Interest
194 cout << "\n>>>> RESULT : " << optMCMC.confLevel * 100 << "% interval on " << firstPOI->GetName() << " is : ["
195 << interval->LowerLimit(*firstPOI) << ", " << interval->UpperLimit(*firstPOI) << "] " << endl;
196
197 gPad = c1;
198}
double sqrt(double)
#define gROOT
Definition: TROOT.h:406
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
#define gPad
Definition: TVirtualPad.h:287
Int_t getSize() const
RooAbsArg * first() const
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
void setMin(const char *name, Double_t value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:466
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:496
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.
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
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
Class implementing a proposal function that samples the parameter space by moving only in one coordin...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
void Print(Option_t *opts=0) const
Print contents of the workspace.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
The Canvas class.
Definition: TCanvas.h:27
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3942
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition: TSystem.cxx:1291
return c1
Definition: legend1.C:41
const Int_t n
Definition: legend1.C:16
return c2
Definition: legend2.C:14
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
Int_t CeilNint(Double_t x)
Definition: TMath.h:689
Definition: file.py:1