ModelInspector.C: RooStats Model Inspector | Roostats tutorials | OneSidedFrequentistUpperLimitWithBands.C: |
// comparison of MCMC and PLC in a multi-variate gaussian problem #include "RooGlobalFunc.h" #include <stdlib.h> #include "TMatrixDSym.h" #include "RooMultiVarGaussian.h" #include "RooArgList.h" #include "RooRealVar.h" #include "TH2F.h" #include "TCanvas.h" #include "RooAbsReal.h" #include "RooFitResult.h" #include "TStopwatch.h" #include "RooStats/MCMCCalculator.h" #include "RooStats/MetropolisHastings.h" #include "RooStats/MarkovChain.h" #include "RooStats/ConfInterval.h" #include "RooStats/MCMCInterval.h" #include "RooStats/MCMCIntervalPlot.h" #include "RooStats/ModelConfig.h" #include "RooStats/ProposalHelper.h" #include "RooStats/ProposalFunction.h" #include "RooStats/PdfProposal.h" #include "RooStats/ProfileLikelihoodCalculator.h" #include "RooStats/LikelihoodIntervalPlot.h" #include "RooStats/LikelihoodInterval.h" using namespace std; using namespace RooFit; using namespace RooStats; void MultivariateGaussianTest(Int_t dim = 4, Int_t nPOI = 2) { /* Authors: Kevin Belasco and Kyle Cranmer. This tutorial produces an N-dimensional multivariate Gaussian with a non-trivial covariance matrix. By default N=4 (called "dim"). A subset of these are considered parameters of interest. This problem is tractable analytically. We use this mainly as a test of Markov Chain Monte Carlo and we compare the result to the profile likelihood ratio. We use the proposal helper to create a customized proposal function for this problem. For N=4 and 2 parameters of interest it takes about 10-20 seconds and the acceptance rate is 37% */ // let's time this challenging example TStopwatch t; t.Start(); RooArgList xVec; RooArgList muVec; RooArgSet poi; // make the observable and means Int_t i,j; RooRealVar* x; RooRealVar* mu_x; for (i = 0; i < dim; i++) { char* name = Form("x%d", i); x = new RooRealVar(name, name, 0, -3,3); xVec.add(*x); char* mu_name = Form("mu_x%d",i); mu_x = new RooRealVar(mu_name, mu_name, 0, -2,2); muVec.add(*mu_x); } // put them into the list of parameters of interest for (i = 0; i < nPOI; i++) { poi.add(*muVec.at(i)); } // make a covariance matrix that is all 1's TMatrixDSym cov(dim); for (i = 0; i < dim; i++) { for (j = 0; j < dim; j++) { if (i == j) cov(i,j) = 3.; else cov(i,j) = 1.0; } } // now make the multivariate Gaussian RooMultiVarGaussian mvg("mvg", "mvg", xVec, muVec, cov); /////////////////////////////////////////// // make a toy dataset RooDataSet* data = mvg.generate(xVec, 100); /////////////////////////////////////////////// // now create the model config for this problem RooWorkspace* w = new RooWorkspace("MVG"); ModelConfig modelConfig(w); modelConfig.SetPdf(mvg); modelConfig.SetParametersOfInterest(poi); /////////////////////////////////////////// // Setup calculators // MCMC // we want to setup an efficient proposal function // using the covariance matrix from a fit to the data RooFitResult* fit = mvg.fitTo(*data, Save(true)); ProposalHelper ph; ph.SetVariables((RooArgSet&)fit->floatParsFinal()); ph.SetCovMatrix(fit->covarianceMatrix()); ph.SetUpdateProposalParameters(true); ph.SetCacheSize(100); ProposalFunction* pdfProp = ph.GetProposalFunction(); // now create the calculator MCMCCalculator mc(*data, modelConfig); mc.SetConfidenceLevel(0.95); mc.SetNumBurnInSteps(100); mc.SetNumIters(10000); mc.SetNumBins(50); mc.SetProposalFunction(*pdfProp); MCMCInterval* mcInt = mc.GetInterval(); RooArgList* poiList = mcInt->GetAxes(); // now setup the profile likelihood calculator ProfileLikelihoodCalculator plc(*data, modelConfig); plc.SetConfidenceLevel(0.95); LikelihoodInterval* plInt = (LikelihoodInterval*)plc.GetInterval(); /////////////////////////////////////////// // make some plots MCMCIntervalPlot mcPlot(*mcInt); TCanvas* c1 = new TCanvas(); mcPlot.SetLineColor(kGreen); mcPlot.SetLineWidth(2); mcPlot.Draw(); LikelihoodIntervalPlot plPlot(plInt); plPlot.Draw("same"); if (poiList->getSize() == 1) { RooRealVar* p = (RooRealVar*)poiList->at(0); Double_t ll = mcInt->LowerLimit(*p); Double_t ul = mcInt->UpperLimit(*p); cout << "MCMC interval: [" << ll << ", " << ul << "]" << endl; } if (poiList->getSize() == 2) { RooRealVar* p0 = (RooRealVar*)poiList->at(0); RooRealVar* p1 = (RooRealVar*)poiList->at(1); TCanvas* scatter = new TCanvas(); Double_t ll = mcInt->LowerLimit(*p0); Double_t ul = mcInt->UpperLimit(*p0); cout << "MCMC interval on p0: [" << ll << ", " << ul << "]" << endl; ll = mcInt->LowerLimit(*p0); ul = mcInt->UpperLimit(*p0); cout << "MCMC interval on p1: [" << ll << ", " << ul << "]" << endl; //MCMC interval on p0: [-0.2, 0.6] //MCMC interval on p1: [-0.2, 0.6] mcPlot.DrawChainScatter(*p0, *p1); scatter->Update(); } t.Print(); } MultivariateGaussianTest.C:1 MultivariateGaussianTest.C:2 MultivariateGaussianTest.C:3 MultivariateGaussianTest.C:4 MultivariateGaussianTest.C:5 MultivariateGaussianTest.C:6 MultivariateGaussianTest.C:7 MultivariateGaussianTest.C:8 MultivariateGaussianTest.C:9 MultivariateGaussianTest.C:10 MultivariateGaussianTest.C:11 MultivariateGaussianTest.C:12 MultivariateGaussianTest.C:13 MultivariateGaussianTest.C:14 MultivariateGaussianTest.C:15 MultivariateGaussianTest.C:16 MultivariateGaussianTest.C:17 MultivariateGaussianTest.C:18 MultivariateGaussianTest.C:19 MultivariateGaussianTest.C:20 MultivariateGaussianTest.C:21 MultivariateGaussianTest.C:22 MultivariateGaussianTest.C:23 MultivariateGaussianTest.C:24 MultivariateGaussianTest.C:25 MultivariateGaussianTest.C:26 MultivariateGaussianTest.C:27 MultivariateGaussianTest.C:28 MultivariateGaussianTest.C:29 MultivariateGaussianTest.C:30 MultivariateGaussianTest.C:31 MultivariateGaussianTest.C:32 MultivariateGaussianTest.C:33 MultivariateGaussianTest.C:34 MultivariateGaussianTest.C:35 MultivariateGaussianTest.C:36 MultivariateGaussianTest.C:37 MultivariateGaussianTest.C:38 MultivariateGaussianTest.C:39 MultivariateGaussianTest.C:40 MultivariateGaussianTest.C:41 MultivariateGaussianTest.C:42 MultivariateGaussianTest.C:43 MultivariateGaussianTest.C:44 MultivariateGaussianTest.C:45 MultivariateGaussianTest.C:46 MultivariateGaussianTest.C:47 MultivariateGaussianTest.C:48 MultivariateGaussianTest.C:49 MultivariateGaussianTest.C:50 MultivariateGaussianTest.C:51 MultivariateGaussianTest.C:52 MultivariateGaussianTest.C:53 MultivariateGaussianTest.C:54 MultivariateGaussianTest.C:55 MultivariateGaussianTest.C:56 MultivariateGaussianTest.C:57 MultivariateGaussianTest.C:58 MultivariateGaussianTest.C:59 MultivariateGaussianTest.C:60 MultivariateGaussianTest.C:61 MultivariateGaussianTest.C:62 MultivariateGaussianTest.C:63 MultivariateGaussianTest.C:64 MultivariateGaussianTest.C:65 MultivariateGaussianTest.C:66 MultivariateGaussianTest.C:67 MultivariateGaussianTest.C:68 MultivariateGaussianTest.C:69 MultivariateGaussianTest.C:70 MultivariateGaussianTest.C:71 MultivariateGaussianTest.C:72 MultivariateGaussianTest.C:73 MultivariateGaussianTest.C:74 MultivariateGaussianTest.C:75 MultivariateGaussianTest.C:76 MultivariateGaussianTest.C:77 MultivariateGaussianTest.C:78 MultivariateGaussianTest.C:79 MultivariateGaussianTest.C:80 MultivariateGaussianTest.C:81 MultivariateGaussianTest.C:82 MultivariateGaussianTest.C:83 MultivariateGaussianTest.C:84 MultivariateGaussianTest.C:85 MultivariateGaussianTest.C:86 MultivariateGaussianTest.C:87 MultivariateGaussianTest.C:88 MultivariateGaussianTest.C:89 MultivariateGaussianTest.C:90 MultivariateGaussianTest.C:91 MultivariateGaussianTest.C:92 MultivariateGaussianTest.C:93 MultivariateGaussianTest.C:94 MultivariateGaussianTest.C:95 MultivariateGaussianTest.C:96 MultivariateGaussianTest.C:97 MultivariateGaussianTest.C:98 MultivariateGaussianTest.C:99 MultivariateGaussianTest.C:100 MultivariateGaussianTest.C:101 MultivariateGaussianTest.C:102 MultivariateGaussianTest.C:103 MultivariateGaussianTest.C:104 MultivariateGaussianTest.C:105 MultivariateGaussianTest.C:106 MultivariateGaussianTest.C:107 MultivariateGaussianTest.C:108 MultivariateGaussianTest.C:109 MultivariateGaussianTest.C:110 MultivariateGaussianTest.C:111 MultivariateGaussianTest.C:112 MultivariateGaussianTest.C:113 MultivariateGaussianTest.C:114 MultivariateGaussianTest.C:115 MultivariateGaussianTest.C:116 MultivariateGaussianTest.C:117 MultivariateGaussianTest.C:118 MultivariateGaussianTest.C:119 MultivariateGaussianTest.C:120 MultivariateGaussianTest.C:121 MultivariateGaussianTest.C:122 MultivariateGaussianTest.C:123 MultivariateGaussianTest.C:124 MultivariateGaussianTest.C:125 MultivariateGaussianTest.C:126 MultivariateGaussianTest.C:127 MultivariateGaussianTest.C:128 MultivariateGaussianTest.C:129 MultivariateGaussianTest.C:130 MultivariateGaussianTest.C:131 MultivariateGaussianTest.C:132 MultivariateGaussianTest.C:133 MultivariateGaussianTest.C:134 MultivariateGaussianTest.C:135 MultivariateGaussianTest.C:136 MultivariateGaussianTest.C:137 MultivariateGaussianTest.C:138 MultivariateGaussianTest.C:139 MultivariateGaussianTest.C:140 MultivariateGaussianTest.C:141 MultivariateGaussianTest.C:142 MultivariateGaussianTest.C:143 MultivariateGaussianTest.C:144 MultivariateGaussianTest.C:145 MultivariateGaussianTest.C:146 MultivariateGaussianTest.C:147 MultivariateGaussianTest.C:148 MultivariateGaussianTest.C:149 MultivariateGaussianTest.C:150 MultivariateGaussianTest.C:151 MultivariateGaussianTest.C:152 MultivariateGaussianTest.C:153 MultivariateGaussianTest.C:154 MultivariateGaussianTest.C:155 MultivariateGaussianTest.C:156 MultivariateGaussianTest.C:157 MultivariateGaussianTest.C:158 MultivariateGaussianTest.C:159 MultivariateGaussianTest.C:160 MultivariateGaussianTest.C:161 MultivariateGaussianTest.C:162 MultivariateGaussianTest.C:163 MultivariateGaussianTest.C:164 MultivariateGaussianTest.C:165 MultivariateGaussianTest.C:166 MultivariateGaussianTest.C:167 MultivariateGaussianTest.C:168 MultivariateGaussianTest.C:169 MultivariateGaussianTest.C:170 MultivariateGaussianTest.C:171 MultivariateGaussianTest.C:172 MultivariateGaussianTest.C:173 MultivariateGaussianTest.C:174 MultivariateGaussianTest.C:175 MultivariateGaussianTest.C:176 |
|