Logo ROOT  
Reference Guide
StandardBayesianNumericalDemo.C File Reference

Detailed Description

Standard demo of the numerical Bayesian calculator.

View in nbviewer Open in SWAN This is a standard demo that can be used with any ROOT file prepared in the standard way. You specify:

  • name for input ROOT file
  • name of workspace inside ROOT file that holds model and data
  • name of ModelConfig that specifies details for calculator tools
  • name of dataset

With default parameters the macro will attempt to run the standard hist2workspace example and read the ROOT file that it produces.

The actual heart of the demo is only about 10 lines long.

The BayesianCalculator is based on Bayes's theorem and performs the integration using ROOT's numeric integration utilities

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:ObjectHandling -- RooWorkspace::import(combined) importing RooUniform::prior
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_simPdf_FOR_OBS_channelCat:obs_x_channel1 with 4 entries
[#1] INFO:Minization -- Including the following constraint terms in minimization: (alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#1] INFO:Minization -- The following global observables have been defined: (nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state channel1 (2 dataset entries)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 1 slave calculators.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Eval -- BayesianCalculator::GetPosteriorFunction : nll value -1033.75 poi value = 1
[#1] INFO:Eval -- BayesianCalculator::GetPosteriorFunction : minimum of NLL vs POI for POI = 1.12121 min NLL = -1033.78
[#1] INFO:InputArguments -- BayesianCalculator:GetInterval Compute the interval from the posterior cdf
[#1] INFO:NumericIntegration -- PosteriorCdfFunction - integral of posterior = 0.0560849 +/- 0.000423427
[#0] WARNING:Eval -- BayesianCalculator::GetInterval : 315 errors reported in evaluating log-likelihood function
[#1] INFO:Eval -- BayesianCalculator::GetInterval - found a valid interval : [0.170773 , 2.33979 ]
>>>> RESULT : 95% interval on SigXsecOverSM is : [0.170773, 2.33979]
#include "TFile.h"
#include "TROOT.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRealVar.h"
#include "RooUniform.h"
#include "RooPlot.h"
#include "TSystem.h"
#include <cassert>
using namespace RooFit;
using namespace RooStats;
struct BayesianNumericalOptions {
double confLevel = 0.95; // interval CL
TString integrationType = ""; // integration Type (default is adaptive (numerical integration)
// possible values are "TOYMC" (toy MC integration, work when nuisances have a constraints pdf)
// "VEGAS" , "MISER", or "PLAIN" (these are all possible MC integration)
int nToys =
10000; // number of toys used for the MC integrations - for Vegas should be probably set to an higher value
bool scanPosterior =
false; // flag to compute interval by scanning posterior (it is more robust but maybe less precise)
bool plotPosterior = false; // plot posterior function after having computed the interval
int nScanPoints = 50; // number of points for scanning the posterior (if scanPosterior = false it is used only for
// plotting). Use by default a low value to speed-up tutorial
int intervalType = 1; // type of interval (0 is shortest, 1 central, 2 upper limit)
double maxPOI = -999; // force a different value of POI for doing the scan (default is given value)
double nSigmaNuisance = -1; // force integration of nuisance parameters to be within nSigma of their error (do first
// a model fit to find nuisance error)
};
BayesianNumericalOptions optBayes;
void StandardBayesianNumericalDemo(const char *infile = "", const char *workspaceName = "combined",
const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
{
// option definitions
double confLevel = optBayes.confLevel;
TString integrationType = optBayes.integrationType;
int nToys = optBayes.nToys;
bool scanPosterior = optBayes.scanPosterior;
bool plotPosterior = optBayes.plotPosterior;
int nScanPoints = optBayes.nScanPoints;
int intervalType = optBayes.intervalType;
int maxPOI = optBayes.maxPOI;
double nSigmaNuisance = optBayes.nSigmaNuisance;
// -------------------------------------------------------
// First part is just to access a user-defined file
// or create the standard example file if it doesn't exist
const char *filename = "";
if (!strcmp(infile, "")) {
filename = "results/example_combined_GaussExample_model.root";
bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
// if file does not exists generate with histfactory
if (!fileExist) {
#ifdef _WIN32
cout << "HistFactory file cannot be generated on Windows - exit" << endl;
return;
#endif
// Normally this would be run on the command line
cout << "will run standard hist2workspace example" << endl;
gROOT->ProcessLine(".! prepareHistFactory .");
gROOT->ProcessLine(".! hist2workspace config/example.xml");
cout << "\n\n---------------------" << endl;
cout << "Done creating example input" << endl;
cout << "---------------------\n\n" << endl;
}
} else
filename = infile;
// Try to open the file
TFile *file = TFile::Open(filename);
// if input file was specified byt not found, quit
if (!file) {
cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return;
}
// -------------------------------------------------------
// Tutorial starts here
// -------------------------------------------------------
// get the workspace out of the file
RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
cout << "workspace not found" << endl;
return;
}
// get the modelConfig out of the file
ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
// get the modelConfig out of the file
RooAbsData *data = w->data(dataName);
// make sure ingredients are found
if (!data || !mc) {
w->Print();
cout << "data or ModelConfig was not found" << endl;
return;
}
// ------------------------------------------
// create and use the BayesianCalculator
// to find and plot the 95% credible interval
// on the parameter of interest as specified
// in the model config
// before we do that, we must specify our prior
// it belongs in the model config, but it may not have
// been specified
RooUniform prior("prior", "", *mc->GetParametersOfInterest());
w->import(prior);
mc->SetPriorPdf(*w->pdf("prior"));
// do without systematics
// mc->SetNuisanceParameters(RooArgSet() );
if (nSigmaNuisance > 0) {
RooAbsPdf *pdf = mc->GetPdf();
assert(pdf);
RooFitResult *res =
res->Print();
for (int i = 0; i < nuisPar.getSize(); ++i) {
RooRealVar *v = dynamic_cast<RooRealVar *>(&nuisPar[i]);
assert(v);
v->setMin(TMath::Max(v->getMin(), v->getVal() - nSigmaNuisance * v->getError()));
v->setMax(TMath::Min(v->getMax(), v->getVal() + nSigmaNuisance * v->getError()));
std::cout << "setting interval for nuisance " << v->GetName() << " : [ " << v->getMin() << " , "
<< v->getMax() << " ]" << std::endl;
}
}
BayesianCalculator bayesianCalc(*data, *mc);
bayesianCalc.SetConfidenceLevel(confLevel); // 95% interval
// default of the calculator is central interval. here use shortest , central or upper limit depending on input
// doing a shortest interval might require a longer time since it requires a scan of the posterior function
if (intervalType == 0)
bayesianCalc.SetShortestInterval(); // for shortest interval
if (intervalType == 1)
bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval
if (intervalType == 2)
bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit
if (!integrationType.IsNull()) {
bayesianCalc.SetIntegrationType(integrationType); // set integrationType
bayesianCalc.SetNumIters(nToys); // set number of iterations (i.e. number of toys for MC integrations)
}
// in case of toyMC make a nuisance pdf
if (integrationType.Contains("TOYMC")) {
RooAbsPdf *nuisPdf = RooStats::MakeNuisancePdf(*mc, "nuisance_pdf");
cout << "using TOYMC integration: make nuisance pdf from the model " << std::endl;
nuisPdf->Print();
bayesianCalc.ForceNuisancePdf(*nuisPdf);
scanPosterior = true; // for ToyMC the posterior is scanned anyway so used given points
}
// compute interval by scanning the posterior function
if (scanPosterior)
bayesianCalc.SetScanOfPosterior(nScanPoints);
if (maxPOI != -999 && maxPOI > poi->getMin())
poi->setMax(maxPOI);
SimpleInterval *interval = bayesianCalc.GetInterval();
// print out the interval on the first Parameter of Interest
cout << "\n>>>> RESULT : " << confLevel * 100 << "% interval on " << poi->GetName() << " is : ["
<< interval->LowerLimit() << ", " << interval->UpperLimit() << "] " << endl;
// end in case plotting is not requested
if (!plotPosterior)
return;
// make a plot
// since plotting may take a long time (it requires evaluating
// the posterior in many points) this command will speed up
// by reducing the number of points to plot - do 50
// ignore errors of PDF if is zero
cout << "\nDrawing plot of posterior function....." << endl;
// always plot using numer of scan points
bayesianCalc.SetScanOfPosterior(nScanPoints);
RooPlot *plot = bayesianCalc.GetPosteriorPlot();
plot->Draw();
}
Author
Kyle Cranmer

Definition in file StandardBayesianNumericalDemo.C.

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
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
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
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:3946
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
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
Definition: RooFitResult.h:40
RooStats::SimpleInterval
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:102
TSystem.h
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
ModelConfig.h
RooFit
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
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
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:246
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:249
RooFit::Minimizer
RooCmdArg Minimizer(const char *type, const char *alg=0)
Definition: RooGlobalFunc.cxx:211
RooStatsUtils.h
RooWorkspace
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:99
RooAbsData.h
RooStats
Definition: Asimov.h:19
file
Definition: file.py:1
ROOT::Math::MinimizerOptions::DefaultPrintLevel
static int DefaultPrintLevel()
Definition: MinimizerOptions.cxx:99
RooStats::BayesianCalculator
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:53
RooStats::ModelConfig::GetNuisanceParameters
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:252
RooRealVar
Definition: RooRealVar.h:36
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:187
RooStats::ModelConfig
Definition: ModelConfig.h:36
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
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