Logo ROOT   6.10/09
Reference Guide
MultivariateGaussianTest.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook
4 /// Comparison of MCMC and PLC in a multi-variate gaussian problem
5 ///
6 /// This tutorial produces an N-dimensional multivariate Gaussian
7 /// with a non-trivial covariance matrix. By default N=4 (called "dim").
8 ///
9 /// A subset of these are considered parameters of interest.
10 /// This problem is tractable analytically.
11 ///
12 /// We use this mainly as a test of Markov Chain Monte Carlo
13 /// and we compare the result to the profile likelihood ratio.
14 ///
15 /// We use the proposal helper to create a customized
16 /// proposal function for this problem.
17 ///
18 /// For N=4 and 2 parameters of interest it takes about 10-20 seconds
19 /// and the acceptance rate is 37%
20 ///
21 /// \macro_image
22 /// \macro_output
23 /// \macro_code
24 ///
25 /// \authors Kevin Belasco and Kyle Cranmer
26 
27 #include "RooGlobalFunc.h"
28 #include <stdlib.h>
29 #include "TMatrixDSym.h"
30 #include "RooMultiVarGaussian.h"
31 #include "RooArgList.h"
32 #include "RooRealVar.h"
33 #include "TH2F.h"
34 #include "TCanvas.h"
35 #include "RooAbsReal.h"
36 #include "RooFitResult.h"
37 #include "TStopwatch.h"
40 #include "RooStats/MarkovChain.h"
41 #include "RooStats/ConfInterval.h"
42 #include "RooStats/MCMCInterval.h"
44 #include "RooStats/ModelConfig.h"
47 #include "RooStats/PdfProposal.h"
51 
52 using namespace std;
53 using namespace RooFit;
54 using namespace RooStats;
55 
56 
57 void MultivariateGaussianTest(Int_t dim = 4, Int_t nPOI = 2)
58 {
59  // let's time this challenging example
60  TStopwatch t;
61  t.Start();
62 
63  RooArgList xVec;
64  RooArgList muVec;
65  RooArgSet poi;
66 
67  // make the observable and means
68  Int_t i,j;
69  RooRealVar* x;
70  RooRealVar* mu_x;
71  for (i = 0; i < dim; i++) {
72  char* name = Form("x%d", i);
73  x = new RooRealVar(name, name, 0, -3,3);
74  xVec.add(*x);
75 
76  char* mu_name = Form("mu_x%d",i);
77  mu_x = new RooRealVar(mu_name, mu_name, 0, -2,2);
78  muVec.add(*mu_x);
79  }
80 
81  // put them into the list of parameters of interest
82  for (i = 0; i < nPOI; i++) {
83  poi.add(*muVec.at(i));
84  }
85 
86  // make a covariance matrix that is all 1's
87  TMatrixDSym cov(dim);
88  for (i = 0; i < dim; i++) {
89  for (j = 0; j < dim; j++) {
90  if (i == j) cov(i,j) = 3.;
91  else cov(i,j) = 1.0;
92  }
93  }
94 
95  // now make the multivariate Gaussian
96  RooMultiVarGaussian mvg("mvg", "mvg", xVec, muVec, cov);
97 
98  // --------------------
99  // make a toy dataset
100  RooDataSet* data = mvg.generate(xVec, 100);
101 
102  // now create the model config for this problem
103  RooWorkspace* w = new RooWorkspace("MVG");
104  ModelConfig modelConfig(w);
105  modelConfig.SetPdf(mvg);
106  modelConfig.SetParametersOfInterest(poi);
107 
108  // -------------------------------------------------------
109  // Setup calculators
110 
111  // MCMC
112  // we want to setup an efficient proposal function
113  // using the covariance matrix from a fit to the data
114  RooFitResult* fit = mvg.fitTo(*data, Save(true));
115  ProposalHelper ph;
117  ph.SetCovMatrix(fit->covarianceMatrix());
119  ph.SetCacheSize(100);
120  ProposalFunction* pdfProp = ph.GetProposalFunction();
121 
122  // now create the calculator
123  MCMCCalculator mc(*data, modelConfig);
124  mc.SetConfidenceLevel(0.95);
125  mc.SetNumBurnInSteps(100);
126  mc.SetNumIters(10000);
127  mc.SetNumBins(50);
128  mc.SetProposalFunction(*pdfProp);
129 
130  MCMCInterval* mcInt = mc.GetInterval();
131  RooArgList* poiList = mcInt->GetAxes();
132 
133  // now setup the profile likelihood calculator
134  ProfileLikelihoodCalculator plc(*data, modelConfig);
135  plc.SetConfidenceLevel(0.95);
136  LikelihoodInterval* plInt = (LikelihoodInterval*)plc.GetInterval();
137 
138 
139  // make some plots
140  MCMCIntervalPlot mcPlot(*mcInt);
141 
142  TCanvas* c1 = new TCanvas();
143  mcPlot.SetLineColor(kGreen);
144  mcPlot.SetLineWidth(2);
145  mcPlot.Draw();
146 
147  LikelihoodIntervalPlot plPlot(plInt);
148  plPlot.Draw("same");
149 
150  if (poiList->getSize() == 1) {
151  RooRealVar* p = (RooRealVar*)poiList->at(0);
152  Double_t ll = mcInt->LowerLimit(*p);
153  Double_t ul = mcInt->UpperLimit(*p);
154  cout << "MCMC interval: [" << ll << ", " << ul << "]" << endl;
155  }
156 
157  if (poiList->getSize() == 2) {
158  RooRealVar* p0 = (RooRealVar*)poiList->at(0);
159  RooRealVar* p1 = (RooRealVar*)poiList->at(1);
160  TCanvas* scatter = new TCanvas();
161  Double_t ll = mcInt->LowerLimit(*p0);
162  Double_t ul = mcInt->UpperLimit(*p0);
163  cout << "MCMC interval on p0: [" << ll << ", " << ul << "]" << endl;
164  ll = mcInt->LowerLimit(*p0);
165  ul = mcInt->UpperLimit(*p0);
166  cout << "MCMC interval on p1: [" << ll << ", " << ul << "]" << endl;
167 
168  //MCMC interval on p0: [-0.2, 0.6]
169  //MCMC interval on p1: [-0.2, 0.6]
170 
171  mcPlot.DrawChainScatter(*p0, *p1);
172  scatter->Update();
173  }
174 
175  t.Print();
176 
177 }
178 
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:108
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:86
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
This class provides simple and straightforward utilities to plot a MCMCInterval object.
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
return c1
Definition: legend1.C:41
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
virtual ProposalFunction * GetProposalFunction()
virtual void SetUpdateProposalParameters(Bool_t updateParams)
This class provides simple and straightforward utilities to plot a LikelihoodInterval object...
ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface class f...
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList ...
Definition: MCMCInterval.h:99
int Int_t
Definition: RtypesCore.h:41
Definition: Rtypes.h:56
STL namespace.
Double_t x[n]
Definition: legend1.C:17
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
Int_t getSize() const
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
virtual void SetCacheSize(Int_t size)
char * Form(const char *fmt,...)
static double p1(double t, double a, double b)
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Multivariate Gaussian p.d.f.
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
double Double_t
Definition: RtypesCore.h:55
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
virtual void SetCovMatrix(const TMatrixDSym &covMatrix)
RooCmdArg Save(Bool_t flag=kTRUE)
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void SetVariables(RooArgList &vars)
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2208
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:42
Stopwatch class.
Definition: TStopwatch.h:28