Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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"
39#include "RooPlot.h"
40#include "TSystem.h"
41
42#include <cassert>
43
44using namespace RooFit;
45using namespace RooStats;
46
47struct 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
66BayesianNumericalOptions optBayes;
67
68void 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 // Normally this would be run on the command line
94 cout << "will run standard hist2workspace example" << endl;
95 gROOT->ProcessLine(".! prepareHistFactory .");
96 gROOT->ProcessLine(".! hist2workspace config/example.xml");
97 cout << "\n\n---------------------" << endl;
98 cout << "Done creating example input" << endl;
99 cout << "---------------------\n\n" << endl;
100 }
101
102 } else
103 filename = infile;
104
105 // Try to open the file
106 TFile *file = TFile::Open(filename);
107
108 // if input file was specified but not found, quit
109 if (!file) {
110 cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
111 return;
112 }
113
114 // -------------------------------------------------------
115 // Tutorial starts here
116 // -------------------------------------------------------
117
118 // get the workspace out of the file
119 RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
120 if (!w) {
121 cout << "workspace not found" << endl;
122 return;
123 }
124
125 // get the modelConfig out of the file
126 ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
127
128 // get the modelConfig out of the file
129 RooAbsData *data = w->data(dataName);
130
131 // make sure ingredients are found
132 if (!data || !mc) {
133 w->Print();
134 cout << "data or ModelConfig was not found" << endl;
135 return;
136 }
137
138 // ------------------------------------------
139 // create and use the BayesianCalculator
140 // to find and plot the 95% credible interval
141 // on the parameter of interest as specified
142 // in the model config
143
144 // before we do that, we must specify our prior
145 // it belongs in the model config, but it may not have
146 // been specified
147 RooUniform prior("prior", "", *mc->GetParametersOfInterest());
148 w->import(prior);
149 mc->SetPriorPdf(*w->pdf("prior"));
150
151 // do without systematics
152 // mc->SetNuisanceParameters(RooArgSet() );
153 if (nSigmaNuisance > 0) {
154 RooAbsPdf *pdf = mc->GetPdf();
155 assert(pdf);
156 std::unique_ptr<RooFitResult> res{
159
160 res->Print();
161 RooArgList nuisPar(*mc->GetNuisanceParameters());
162 for (int i = 0; i < nuisPar.getSize(); ++i) {
163 RooRealVar *v = dynamic_cast<RooRealVar *>(&nuisPar[i]);
164 assert(v);
165 v->setMin(TMath::Max(v->getMin(), v->getVal() - nSigmaNuisance * v->getError()));
166 v->setMax(TMath::Min(v->getMax(), v->getVal() + nSigmaNuisance * v->getError()));
167 std::cout << "setting interval for nuisance " << v->GetName() << " : [ " << v->getMin() << " , "
168 << v->getMax() << " ]" << std::endl;
169 }
170 }
171
172 BayesianCalculator bayesianCalc(*data, *mc);
173 bayesianCalc.SetConfidenceLevel(confLevel); // 95% interval
174
175 // default of the calculator is central interval. here use shortest , central or upper limit depending on input
176 // doing a shortest interval might require a longer time since it requires a scan of the posterior function
177 if (intervalType == 0)
178 bayesianCalc.SetShortestInterval(); // for shortest interval
179 if (intervalType == 1)
180 bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval
181 if (intervalType == 2)
182 bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit
183
184 if (!integrationType.IsNull()) {
185 bayesianCalc.SetIntegrationType(integrationType); // set integrationType
186 bayesianCalc.SetNumIters(nToys); // set number of iterations (i.e. number of toys for MC integrations)
187 }
188
189 // in case of toyMC make a nuisance pdf
190 if (integrationType.Contains("TOYMC")) {
191 RooAbsPdf *nuisPdf = RooStats::MakeNuisancePdf(*mc, "nuisance_pdf");
192 cout << "using TOYMC integration: make nuisance pdf from the model " << std::endl;
193 nuisPdf->Print();
194 bayesianCalc.ForceNuisancePdf(*nuisPdf);
195 scanPosterior = true; // for ToyMC the posterior is scanned anyway so used given points
196 }
197
198 // compute interval by scanning the posterior function
199 if (scanPosterior)
200 bayesianCalc.SetScanOfPosterior(nScanPoints);
201
203 if (maxPOI != -999 && maxPOI > poi->getMin())
204 poi->setMax(maxPOI);
205
206 SimpleInterval *interval = bayesianCalc.GetInterval();
207
208 // print out the interval on the first Parameter of Interest
209 cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << poi->GetName() << " is : ["
210 << interval->LowerLimit() << ", " << interval->UpperLimit() << "] " << endl;
211
212 // end in case plotting is not requested
213 if (!plotPosterior)
214 return;
215
216 // make a plot
217 // since plotting may take a long time (it requires evaluating
218 // the posterior in many points) this command will speed up
219 // by reducing the number of points to plot - do 50
220
221 // ignore errors of PDF if is zero
223
224 cout << "\nDrawing plot of posterior function....." << endl;
225
226 // always plot using number of scan points
227 bayesianCalc.SetScanOfPosterior(nScanPoints);
228
229 RooPlot *plot = bayesianCalc.GetPosteriorPlot();
230 plot->Draw();
231}
winID h TVirtualViewer3D TVirtualGLPainter char TVirtualGLPainter plot
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
#define gROOT
Definition TROOT.h:406
R__EXTERN TSystem * gSystem
Definition TSystem.h:561
static const std::string & DefaultMinimizerType()
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:263
RooAbsArg * first() const
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:157
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
static void setEvalErrorLoggingMode(ErrorLoggingMode m)
Set evaluation error logging mode.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
Plot frame and a container for graphics objects within that frame.
Definition RooPlot.h:43
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setMax(const char *name, double value)
Set maximum of name range to given value.
BayesianCalculator is a concrete implementation of IntervalCalculator, providing the computation of a...
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:35
virtual void SetPriorPdf(const RooAbsPdf &pdf)
Set the Prior Pdf, add to the workspace if not already there.
Definition ModelConfig.h:95
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return nullptr if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return nullptr if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return nullptr if pdf has not been specified or does not exist)
SimpleInterval is a concrete implementation of the ConfInterval interface.
virtual double UpperLimit()
return the interval upper limit
virtual double LowerLimit()
return the interval lower limit
Flat p.d.f.
Definition RooUniform.h:24
Persistable container for RooFit projects.
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
A ROOT file is an on-disk file, usually with extension .root, that stores objects in a file-system-li...
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:4086
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:292
Basic string class.
Definition TString.h:139
Bool_t IsNull() const
Definition TString.h:414
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
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:1296
RooCmdArg Minimizer(const char *type, const char *alg=nullptr)
RooCmdArg Hesse(bool flag=true)
RooCmdArg Save(bool flag=true)
RooCmdArg PrintLevel(Int_t code)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:64
Namespace for the RooStats classes.
Definition CodegenImpl.h:58
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
extract constraint terms from pdf
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198