Logo ROOT   6.21/01
Reference Guide
StandardBayesianNumericalDemo.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook -js
4 /// 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 
112  // if input file was specified byt not found, quit
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) {
125  cout << "workspace not found" << endl;
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();
138  cout << "data or ModelConfig was not found" << endl;
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 }
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
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:1279
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
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:3925
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
RooCmdArg PrintLevel(Int_t code)
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:48
#define gROOT
Definition: TROOT.h:415
Basic string class.
Definition: TString.h:131
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:465
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:66
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:281
virtual Double_t LowerLimit()
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
R__EXTERN TSystem * gSystem
Definition: TSystem.h:557
static const std::string & DefaultMinimizerType()
RooAbsArg * first() const
RooCmdArg Minimizer(const char *type, const char *alg=0)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
Flat p.d.f.
Definition: RooUniform.h:24
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
RooCmdArg Hesse(Bool_t flag=kTRUE)
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:87
Bool_t IsNull() const
Definition: TString.h:402
RooCmdArg Save(Bool_t flag=kTRUE)
virtual Double_t UpperLimit()
SimpleInterval is a concrete implementation of the ConfInterval interface.
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
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: file.py:1
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
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:1256
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
void Print(Option_t *opts=0) const
Print contents of the workspace.
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712