From $ROOTSYS/tutorials/roostats/StandardBayesianNumericalDemo.C

// 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"
#include "TSystem.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 = 20;          // 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 withing nSigma of their error (do first a model fit to find nuisance error)

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
  ////////////////////////////////////////////////////////////

   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 = pdf->fitTo(*data, Save(true), Minimizer(ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str()), Hesse(true),
                                     PrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()-1) );

     res->Print();
     RooArgList nuisPar(*mc->GetNuisanceParameters());
     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(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
 StandardBayesianNumericalDemo.C:192
 StandardBayesianNumericalDemo.C:193
 StandardBayesianNumericalDemo.C:194
 StandardBayesianNumericalDemo.C:195
 StandardBayesianNumericalDemo.C:196
 StandardBayesianNumericalDemo.C:197
 StandardBayesianNumericalDemo.C:198
 StandardBayesianNumericalDemo.C:199
 StandardBayesianNumericalDemo.C:200
 StandardBayesianNumericalDemo.C:201
 StandardBayesianNumericalDemo.C:202
 StandardBayesianNumericalDemo.C:203
 StandardBayesianNumericalDemo.C:204
 StandardBayesianNumericalDemo.C:205
 StandardBayesianNumericalDemo.C:206
 StandardBayesianNumericalDemo.C:207