StandardBayesianMCMCDemo.C: Standard demo of the Bayesian MCMC calculator | Roostats tutorials | StandardFeldmanCousinsDemo.C: Standard demo of the Feldman-Cousins calculator |
// Standard demo of the numerical Bayesian calculator /* Author: Kyle Cranmer date: Dec. 2010 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 */ #include "TFile.h" #include "TROOT.h" #include "RooWorkspace.h" #include "RooAbsData.h" #include "RooRealVar.h" #include "RooUniform.h" #include "RooStats/ModelConfig.h" #include "RooStats/BayesianCalculator.h" #include "RooStats/SimpleInterval.h" #include "RooStats/RooStatsUtils.h" #include "RooPlot.h" using namespace RooFit; using namespace RooStats; 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) int nScanPoints = 50; // number of points for scanning the posterior (if scanPosterior = false it is used only for plotting) 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) void StandardBayesianNumericalDemo(const char* infile = "", const char* workspaceName = "combined", const char* modelConfigName = "ModelConfig", const char* dataName = "obsData"){ ///////////////////////////////////////////////////////////// // First part is just to access a user-defined file // or create the standard example file if it doesn't exist //////////////////////////////////////////////////////////// TString filename = infile; if (filename.IsNull()) { filename = "results/example_combined_GaussExample_model.root"; std::cout << "Use standard file generated with HistFactory : " << filename << std::endl; } // Check if example input file exists TFile *file = TFile::Open(filename); // if input file was specified but not found, quit if(!file && !TString(infile).IsNull()){ cout <<"file " << filename << " not found" << endl; return; } // if default file not found, try to create it if(!file ){ // 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; // now try to access the file again file = TFile::Open(filename); } if(!file){ // if it is still not there, then we can't continue cout << "Not able to run hist2workspace to create example input" <<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() ); BayesianCalculator bayesianCalc(*data,*mc); bayesianCalc.SetConfidenceLevel(0.95); // 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 ietrations (i.e. number of toys for MC integrations) } // in case of toyMC make a nnuisance 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); RooRealVar* poi = (RooRealVar*) mc->GetParametersOfInterest()->first(); if (maxPOI != -999 && maxPOI > poi->getMin()) poi->setMax(maxPOI); SimpleInterval* interval = bayesianCalc.GetInterval(); // print out the iterval on the first Parameter of Interest cout << "\n95% interval on " << poi->GetName()<<" is : ["<< interval->LowerLimit() << ", "<< interval->UpperLimit() <<"] "<<endl; // 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 cout << "\nDrawing plot of posterior function....." << endl; bayesianCalc.SetScanOfPosterior(nScanPoints); RooPlot * plot = bayesianCalc.GetPosteriorPlot(); plot->Draw(); } StandardBayesianNumericalDemo.C:1 StandardBayesianNumericalDemo.C:2 StandardBayesianNumericalDemo.C:3 StandardBayesianNumericalDemo.C:4 StandardBayesianNumericalDemo.C:5 StandardBayesianNumericalDemo.C:6 StandardBayesianNumericalDemo.C:7 StandardBayesianNumericalDemo.C:8 StandardBayesianNumericalDemo.C:9 StandardBayesianNumericalDemo.C:10 StandardBayesianNumericalDemo.C:11 StandardBayesianNumericalDemo.C:12 StandardBayesianNumericalDemo.C:13 StandardBayesianNumericalDemo.C:14 StandardBayesianNumericalDemo.C:15 StandardBayesianNumericalDemo.C:16 StandardBayesianNumericalDemo.C:17 StandardBayesianNumericalDemo.C:18 StandardBayesianNumericalDemo.C:19 StandardBayesianNumericalDemo.C:20 StandardBayesianNumericalDemo.C:21 StandardBayesianNumericalDemo.C:22 StandardBayesianNumericalDemo.C:23 StandardBayesianNumericalDemo.C:24 StandardBayesianNumericalDemo.C:25 StandardBayesianNumericalDemo.C:26 StandardBayesianNumericalDemo.C:27 StandardBayesianNumericalDemo.C:28 StandardBayesianNumericalDemo.C:29 StandardBayesianNumericalDemo.C:30 StandardBayesianNumericalDemo.C:31 StandardBayesianNumericalDemo.C:32 StandardBayesianNumericalDemo.C:33 StandardBayesianNumericalDemo.C:34 StandardBayesianNumericalDemo.C:35 StandardBayesianNumericalDemo.C:36 StandardBayesianNumericalDemo.C:37 StandardBayesianNumericalDemo.C:38 StandardBayesianNumericalDemo.C:39 StandardBayesianNumericalDemo.C:40 StandardBayesianNumericalDemo.C:41 StandardBayesianNumericalDemo.C:42 StandardBayesianNumericalDemo.C:43 StandardBayesianNumericalDemo.C:44 StandardBayesianNumericalDemo.C:45 StandardBayesianNumericalDemo.C:46 StandardBayesianNumericalDemo.C:47 StandardBayesianNumericalDemo.C:48 StandardBayesianNumericalDemo.C:49 StandardBayesianNumericalDemo.C:50 StandardBayesianNumericalDemo.C:51 StandardBayesianNumericalDemo.C:52 StandardBayesianNumericalDemo.C:53 StandardBayesianNumericalDemo.C:54 StandardBayesianNumericalDemo.C:55 StandardBayesianNumericalDemo.C:56 StandardBayesianNumericalDemo.C:57 StandardBayesianNumericalDemo.C:58 StandardBayesianNumericalDemo.C:59 StandardBayesianNumericalDemo.C:60 StandardBayesianNumericalDemo.C:61 StandardBayesianNumericalDemo.C:62 StandardBayesianNumericalDemo.C:63 StandardBayesianNumericalDemo.C:64 StandardBayesianNumericalDemo.C:65 StandardBayesianNumericalDemo.C:66 StandardBayesianNumericalDemo.C:67 StandardBayesianNumericalDemo.C:68 StandardBayesianNumericalDemo.C:69 StandardBayesianNumericalDemo.C:70 StandardBayesianNumericalDemo.C:71 StandardBayesianNumericalDemo.C:72 StandardBayesianNumericalDemo.C:73 StandardBayesianNumericalDemo.C:74 StandardBayesianNumericalDemo.C:75 StandardBayesianNumericalDemo.C:76 StandardBayesianNumericalDemo.C:77 StandardBayesianNumericalDemo.C:78 StandardBayesianNumericalDemo.C:79 StandardBayesianNumericalDemo.C:80 StandardBayesianNumericalDemo.C:81 StandardBayesianNumericalDemo.C:82 StandardBayesianNumericalDemo.C:83 StandardBayesianNumericalDemo.C:84 StandardBayesianNumericalDemo.C:85 StandardBayesianNumericalDemo.C:86 StandardBayesianNumericalDemo.C:87 StandardBayesianNumericalDemo.C:88 StandardBayesianNumericalDemo.C:89 StandardBayesianNumericalDemo.C:90 StandardBayesianNumericalDemo.C:91 StandardBayesianNumericalDemo.C:92 StandardBayesianNumericalDemo.C:93 StandardBayesianNumericalDemo.C:94 StandardBayesianNumericalDemo.C:95 StandardBayesianNumericalDemo.C:96 StandardBayesianNumericalDemo.C:97 StandardBayesianNumericalDemo.C:98 StandardBayesianNumericalDemo.C:99 StandardBayesianNumericalDemo.C:100 StandardBayesianNumericalDemo.C:101 StandardBayesianNumericalDemo.C:102 StandardBayesianNumericalDemo.C:103 StandardBayesianNumericalDemo.C:104 StandardBayesianNumericalDemo.C:105 StandardBayesianNumericalDemo.C:106 StandardBayesianNumericalDemo.C:107 StandardBayesianNumericalDemo.C:108 StandardBayesianNumericalDemo.C:109 StandardBayesianNumericalDemo.C:110 StandardBayesianNumericalDemo.C:111 StandardBayesianNumericalDemo.C:112 StandardBayesianNumericalDemo.C:113 StandardBayesianNumericalDemo.C:114 StandardBayesianNumericalDemo.C:115 StandardBayesianNumericalDemo.C:116 StandardBayesianNumericalDemo.C:117 StandardBayesianNumericalDemo.C:118 StandardBayesianNumericalDemo.C:119 StandardBayesianNumericalDemo.C:120 StandardBayesianNumericalDemo.C:121 StandardBayesianNumericalDemo.C:122 StandardBayesianNumericalDemo.C:123 StandardBayesianNumericalDemo.C:124 StandardBayesianNumericalDemo.C:125 StandardBayesianNumericalDemo.C:126 StandardBayesianNumericalDemo.C:127 StandardBayesianNumericalDemo.C:128 StandardBayesianNumericalDemo.C:129 StandardBayesianNumericalDemo.C:130 StandardBayesianNumericalDemo.C:131 StandardBayesianNumericalDemo.C:132 StandardBayesianNumericalDemo.C:133 StandardBayesianNumericalDemo.C:134 StandardBayesianNumericalDemo.C:135 StandardBayesianNumericalDemo.C:136 StandardBayesianNumericalDemo.C:137 StandardBayesianNumericalDemo.C:138 StandardBayesianNumericalDemo.C:139 StandardBayesianNumericalDemo.C:140 StandardBayesianNumericalDemo.C:141 StandardBayesianNumericalDemo.C:142 StandardBayesianNumericalDemo.C:143 StandardBayesianNumericalDemo.C:144 StandardBayesianNumericalDemo.C:145 StandardBayesianNumericalDemo.C:146 StandardBayesianNumericalDemo.C:147 StandardBayesianNumericalDemo.C:148 StandardBayesianNumericalDemo.C:149 StandardBayesianNumericalDemo.C:150 StandardBayesianNumericalDemo.C:151 StandardBayesianNumericalDemo.C:152 StandardBayesianNumericalDemo.C:153 StandardBayesianNumericalDemo.C:154 StandardBayesianNumericalDemo.C:155 StandardBayesianNumericalDemo.C:156 StandardBayesianNumericalDemo.C:157 StandardBayesianNumericalDemo.C:158 StandardBayesianNumericalDemo.C:159 StandardBayesianNumericalDemo.C:160 StandardBayesianNumericalDemo.C:161 StandardBayesianNumericalDemo.C:162 StandardBayesianNumericalDemo.C:163 StandardBayesianNumericalDemo.C:164 StandardBayesianNumericalDemo.C:165 StandardBayesianNumericalDemo.C:166 StandardBayesianNumericalDemo.C:167 StandardBayesianNumericalDemo.C:168 StandardBayesianNumericalDemo.C:169 StandardBayesianNumericalDemo.C:170 StandardBayesianNumericalDemo.C:171 StandardBayesianNumericalDemo.C:172 StandardBayesianNumericalDemo.C:173 StandardBayesianNumericalDemo.C:174 StandardBayesianNumericalDemo.C:175 StandardBayesianNumericalDemo.C:176 StandardBayesianNumericalDemo.C:177 StandardBayesianNumericalDemo.C:178 StandardBayesianNumericalDemo.C:179 StandardBayesianNumericalDemo.C:180 StandardBayesianNumericalDemo.C:181 StandardBayesianNumericalDemo.C:182 StandardBayesianNumericalDemo.C:183 StandardBayesianNumericalDemo.C:184 StandardBayesianNumericalDemo.C:185 StandardBayesianNumericalDemo.C:186 StandardBayesianNumericalDemo.C:187 StandardBayesianNumericalDemo.C:188 StandardBayesianNumericalDemo.C:189 StandardBayesianNumericalDemo.C:190 StandardBayesianNumericalDemo.C:191 |
|