ROOT   6.21/01 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 void MultivariateGaussianTest(Int_t dim = 4, Int_t nPOI = 2)
57 {
58  // let's time this challenging example
59  TStopwatch t;
60  t.Start();
61
62  RooArgList xVec;
63  RooArgList muVec;
64  RooArgSet poi;
65
66  // make the observable and means
67  Int_t i, j;
68  RooRealVar *x;
69  RooRealVar *mu_x;
70  for (i = 0; i < dim; i++) {
71  char *name = Form("x%d", i);
72  x = new RooRealVar(name, name, 0, -3, 3);
74
75  char *mu_name = Form("mu_x%d", i);
76  mu_x = new RooRealVar(mu_name, mu_name, 0, -2, 2);
78  }
79
80  // put them into the list of parameters of interest
81  for (i = 0; i < nPOI; i++) {
83  }
84
85  // make a covariance matrix that is all 1's
86  TMatrixDSym cov(dim);
87  for (i = 0; i < dim; i++) {
88  for (j = 0; j < dim; j++) {
89  if (i == j)
90  cov(i, j) = 3.;
91  else
92  cov(i, j) = 1.0;
93  }
94  }
95
96  // now make the multivariate Gaussian
97  RooMultiVarGaussian mvg("mvg", "mvg", xVec, muVec, cov);
98
99  // --------------------
100  // make a toy dataset
101  RooDataSet *data = mvg.generate(xVec, 100);
102
103  // now create the model config for this problem
104  RooWorkspace *w = new RooWorkspace("MVG");
105  ModelConfig modelConfig(w);
106  modelConfig.SetPdf(mvg);
107  modelConfig.SetParametersOfInterest(poi);
108
109  // -------------------------------------------------------
110  // Setup calculators
111
112  // MCMC
113  // we want to setup an efficient proposal function
114  // using the covariance matrix from a fit to the data
115  RooFitResult *fit = mvg.fitTo(*data, Save(true));
116  ProposalHelper ph;
117  ph.SetVariables((RooArgSet &)fit->floatParsFinal());
118  ph.SetCovMatrix(fit->covarianceMatrix());
120  ph.SetCacheSize(100);
121  ProposalFunction *pdfProp = ph.GetProposalFunction();
122
123  // now create the calculator
124  MCMCCalculator mc(*data, modelConfig);
125  mc.SetConfidenceLevel(0.95);
126  mc.SetNumBurnInSteps(100);
127  mc.SetNumIters(10000);
128  mc.SetNumBins(50);
129  mc.SetProposalFunction(*pdfProp);
130
131  MCMCInterval *mcInt = mc.GetInterval();
132  RooArgList *poiList = mcInt->GetAxes();
133
134  // now setup the profile likelihood calculator
135  ProfileLikelihoodCalculator plc(*data, modelConfig);
136  plc.SetConfidenceLevel(0.95);
137  LikelihoodInterval *plInt = (LikelihoodInterval *)plc.GetInterval();
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 }
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:110
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:88
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
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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...
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList ...
Definition: MCMCInterval.h:99
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
int Int_t
Definition: RtypesCore.h:41
Definition: Rtypes.h:64
STL namespace.
Double_t x[n]
Definition: legend1.C:17
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
RooRealVar represents a fundamental (non-derived) real-valued object.
Definition: RooRealVar.h:36
Int_t getSize() const
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
virtual void SetCacheSize(Int_t size)
char * Form(const char *fmt,...)
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:31
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)
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:30
virtual void SetVariables(RooArgList &vars)
virtual void Update()