ROOT   Reference Guide
StandardBayesianNumericalDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook -js
4 /// \brief Standard demo of the numerical Bayesian 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 BayesianCalculator is based on Bayes's theorem
20 /// and performs the integration using ROOT's numeric integration utilities
21 ///
22 /// \macro_image
23 /// \macro_output
24 /// \macro_code
25 ///
26 /// \author Kyle Cranmer
27
28 #include "TFile.h"
29 #include "TROOT.h"
30 #include "RooWorkspace.h"
31 #include "RooAbsData.h"
32 #include "RooRealVar.h"
33
34 #include "RooUniform.h"
35 #include "RooStats/ModelConfig.h"
38 #include "RooStats/RooStatsUtils.h"
39 #include "RooPlot.h"
40 #include "TSystem.h"
41
42 #include <cassert>
43
44 using namespace RooFit;
45 using namespace RooStats;
46
47 struct BayesianNumericalOptions {
48
49  double confLevel = 0.95; // interval CL
50  TString integrationType = ""; // integration Type (default is adaptive (numerical integration)
51  // possible values are "TOYMC" (toy MC integration, work when nuisances have a constraints pdf)
52  // "VEGAS" , "MISER", or "PLAIN" (these are all possible MC integration)
53  int nToys =
54  10000; // number of toys used for the MC integrations - for Vegas should be probably set to an higher value
55  bool scanPosterior =
56  false; // flag to compute interval by scanning posterior (it is more robust but maybe less precise)
57  bool plotPosterior = false; // plot posterior function after having computed the interval
58  int nScanPoints = 50; // number of points for scanning the posterior (if scanPosterior = false it is used only for
59  // plotting). Use by default a low value to speed-up tutorial
60  int intervalType = 1; // type of interval (0 is shortest, 1 central, 2 upper limit)
61  double maxPOI = -999; // force a different value of POI for doing the scan (default is given value)
62  double nSigmaNuisance = -1; // force integration of nuisance parameters to be within nSigma of their error (do first
63  // a model fit to find nuisance error)
64 };
65
66 BayesianNumericalOptions optBayes;
67
68 void StandardBayesianNumericalDemo(const char *infile = "", const char *workspaceName = "combined",
69  const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
70 {
71
72  // option definitions
73  double confLevel = optBayes.confLevel;
74  TString integrationType = optBayes.integrationType;
75  int nToys = optBayes.nToys;
76  bool scanPosterior = optBayes.scanPosterior;
77  bool plotPosterior = optBayes.plotPosterior;
78  int nScanPoints = optBayes.nScanPoints;
79  int intervalType = optBayes.intervalType;
80  int maxPOI = optBayes.maxPOI;
81  double nSigmaNuisance = optBayes.nSigmaNuisance;
82
83  // -------------------------------------------------------
84  // First part is just to access a user-defined file
85  // or create the standard example file if it doesn't exist
86
87  const char *filename = "";
88  if (!strcmp(infile, "")) {
89  filename = "results/example_combined_GaussExample_model.root";
90  bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
91  // if file does not exists generate with histfactory
92  if (!fileExist) {
93 #ifdef _WIN32
94  cout << "HistFactory file cannot be generated on Windows - exit" << endl;
95  return;
96 #endif
97  // Normally this would be run on the command line
98  cout << "will run standard hist2workspace example" << endl;
99  gROOT->ProcessLine(".! prepareHistFactory .");
100  gROOT->ProcessLine(".! hist2workspace config/example.xml");
101  cout << "\n\n---------------------" << endl;
102  cout << "Done creating example input" << endl;
103  cout << "---------------------\n\n" << endl;
104  }
105
106  } else
107  filename = infile;
108
109  // Try to open the file
110  TFile *file = TFile::Open(filename);
111
113  if (!file) {
114  cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
115  return;
116  }
117
118  // -------------------------------------------------------
119  // Tutorial starts here
120  // -------------------------------------------------------
121
122  // get the workspace out of the file
123  RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
124  if (!w) {
126  return;
127  }
128
129  // get the modelConfig out of the file
130  ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
131
132  // get the modelConfig out of the file
133  RooAbsData *data = w->data(dataName);
134
135  // make sure ingredients are found
136  if (!data || !mc) {
137  w->Print();
139  return;
140  }
141
142  // ------------------------------------------
143  // create and use the BayesianCalculator
144  // to find and plot the 95% credible interval
145  // on the parameter of interest as specified
146  // in the model config
147
148  // before we do that, we must specify our prior
149  // it belongs in the model config, but it may not have
150  // been specified
151  RooUniform prior("prior", "", *mc->GetParametersOfInterest());
152  w->import(prior);
153  mc->SetPriorPdf(*w->pdf("prior"));
154
155  // do without systematics
156  // mc->SetNuisanceParameters(RooArgSet() );
157  if (nSigmaNuisance > 0) {
158  RooAbsPdf *pdf = mc->GetPdf();
159  assert(pdf);
160  RooFitResult *res =
163
164  res->Print();
165  RooArgList nuisPar(*mc->GetNuisanceParameters());
166  for (int i = 0; i < nuisPar.getSize(); ++i) {
167  RooRealVar *v = dynamic_cast<RooRealVar *>(&nuisPar[i]);
168  assert(v);
169  v->setMin(TMath::Max(v->getMin(), v->getVal() - nSigmaNuisance * v->getError()));
170  v->setMax(TMath::Min(v->getMax(), v->getVal() + nSigmaNuisance * v->getError()));
171  std::cout << "setting interval for nuisance " << v->GetName() << " : [ " << v->getMin() << " , "
172  << v->getMax() << " ]" << std::endl;
173  }
174  }
175
176  BayesianCalculator bayesianCalc(*data, *mc);
177  bayesianCalc.SetConfidenceLevel(confLevel); // 95% interval
178
179  // default of the calculator is central interval. here use shortest , central or upper limit depending on input
180  // doing a shortest interval might require a longer time since it requires a scan of the posterior function
181  if (intervalType == 0)
182  bayesianCalc.SetShortestInterval(); // for shortest interval
183  if (intervalType == 1)
184  bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval
185  if (intervalType == 2)
186  bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit
187
188  if (!integrationType.IsNull()) {
189  bayesianCalc.SetIntegrationType(integrationType); // set integrationType
190  bayesianCalc.SetNumIters(nToys); // set number of iterations (i.e. number of toys for MC integrations)
191  }
192
193  // in case of toyMC make a nuisance pdf
194  if (integrationType.Contains("TOYMC")) {
195  RooAbsPdf *nuisPdf = RooStats::MakeNuisancePdf(*mc, "nuisance_pdf");
196  cout << "using TOYMC integration: make nuisance pdf from the model " << std::endl;
197  nuisPdf->Print();
198  bayesianCalc.ForceNuisancePdf(*nuisPdf);
199  scanPosterior = true; // for ToyMC the posterior is scanned anyway so used given points
200  }
201
202  // compute interval by scanning the posterior function
203  if (scanPosterior)
204  bayesianCalc.SetScanOfPosterior(nScanPoints);
205
207  if (maxPOI != -999 && maxPOI > poi->getMin())
208  poi->setMax(maxPOI);
209
210  SimpleInterval *interval = bayesianCalc.GetInterval();
211
212  // print out the interval on the first Parameter of Interest
213  cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << poi->GetName() << " is : ["
214  << interval->LowerLimit() << ", " << interval->UpperLimit() << "] " << endl;
215
216  // end in case plotting is not requested
217  if (!plotPosterior)
218  return;
219
220  // make a plot
221  // since plotting may take a long time (it requires evaluating
222  // the posterior in many points) this command will speed up
223  // by reducing the number of points to plot - do 50
224
225  // ignore errors of PDF if is zero
227
228  cout << "\nDrawing plot of posterior function....." << endl;
229
230  // always plot using numer of scan points
231  bayesianCalc.SetScanOfPosterior(nScanPoints);
232
233  RooPlot *plot = bayesianCalc.GetPosteriorPlot();
234  plot->Draw();
235 }
RooWorkspace::data
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1368
RooUniform
Flat p.d.f.
Definition: RooUniform.h:24
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
RooAbsReal::setEvalErrorLoggingMode
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
Definition: RooAbsReal.cxx:4838
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooAbsData
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:46
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition: RooGlobalFunc.cxx:189
RooAbsArg::Print
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:320
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooFit::Hesse
RooCmdArg Hesse(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:193
TFile::Open
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:3995
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
RooStats::SimpleInterval::LowerLimit
virtual Double_t LowerLimit()
Definition: SimpleInterval.h:47
TString
Basic string class.
Definition: TString.h:136
TSystem::AccessPathName
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:1294
v
@ v
Definition: rootcling_impl.cxx:3635
TFile.h
RooStats::SimpleInterval::UpperLimit
virtual Double_t UpperLimit()
Definition: SimpleInterval.h:49
RooAbsRealLValue::getMin
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
Definition: RooAbsRealLValue.h:86
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooStats::SimpleInterval
SimpleInterval is a concrete implementation of the ConfInterval interface.
Definition: SimpleInterval.h:20
RooWorkspace::import
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
Definition: RooWorkspace.cxx:361
TROOT.h
ROOT::Math::MinimizerOptions::DefaultMinimizerType
static const std::string & DefaultMinimizerType()
Definition: MinimizerOptions.cxx:92
TSystem.h
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
ModelConfig.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
BayesianCalculator.h
RooPlot.h
SimpleInterval.h
RooAbsReal::Ignore
@ Ignore
Definition: RooAbsReal.h:298
RooWorkspace::obj
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
Definition: RooWorkspace.cxx:2106
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
RooWorkspace::pdf
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1277
RooRealVar.h
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
RooStats::ModelConfig::GetPdf
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
RooStats::ModelConfig::GetParametersOfInterest
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
RooFit::Minimizer
RooCmdArg Minimizer(const char *type, const char *alg=0)
Definition: RooGlobalFunc.cxx:211
RooStatsUtils.h
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooRealVar::setMax
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:450
TString::IsNull
Bool_t IsNull() const
Definition: TString.h:407
RooStats::ModelConfig::SetPriorPdf
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:87
RooAbsData.h
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
file
Definition: file.py:1
ROOT::Math::MinimizerOptions::DefaultPrintLevel
static int DefaultPrintLevel()
Definition: MinimizerOptions.cxx:89
RooStats::BayesianCalculator
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
Definition: BayesianCalculator.h:37
RooFitResult::Print
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
RooAbsPdf
Definition: RooAbsPdf.h:40
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooStats::ModelConfig::GetNuisanceParameters
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:36
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
RooUniform.h
RooStats::MakeNuisancePdf
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
Definition: RooStatsUtils.cxx:118
gROOT
#define gROOT
Definition: TROOT.h:406
RooAbsPdf::fitTo
virtual RooFitResult * fitTo(RooAbsData &data, 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())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1261