ROOT   Reference Guide
MultivariateGaussianTest.C File Reference

## Detailed Description

Comparison of MCMC and PLC in a multi-variate gaussian problem

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%

[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
**********
** 1 **SET PRINT 1
**********
**********
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 mu_x0 0.00000e+00 4.00000e-01 -2.00000e+00 2.00000e+00
2 mu_x1 0.00000e+00 4.00000e-01 -2.00000e+00 2.00000e+00
3 mu_x2 0.00000e+00 4.00000e-01 -2.00000e+00 2.00000e+00
4 mu_x3 0.00000e+00 4.00000e-01 -2.00000e+00 2.00000e+00
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 1
**********
**********
** 5 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
FCN=707.822 FROM MIGRAD STATUS=INITIATE 14 CALLS 15 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 mu_x0 0.00000e+00 4.00000e-01 2.01358e-01 -9.83942e+00
2 mu_x1 0.00000e+00 4.00000e-01 2.01358e-01 -1.25127e+01
3 mu_x2 0.00000e+00 4.00000e-01 2.01358e-01 9.88572e+00
4 mu_x3 0.00000e+00 4.00000e-01 2.01358e-01 -4.09155e+00
ERR DEF= 0.5
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=706.56 FROM MIGRAD STATUS=CONVERGED 65 CALLS 66 TOTAL
EDM=1.95881e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 mu_x0 1.80728e-01 1.72980e-01 1.42731e-03 -2.59592e-02
2 mu_x1 2.07351e-01 1.72978e-01 1.42884e-03 -3.69006e-02
3 mu_x2 -1.59412e-02 1.72984e-01 1.42230e-03 3.21539e-02
4 mu_x3 1.23430e-01 1.72982e-01 1.42472e-03 -8.04153e-03
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 4 ERR DEF=0.5
3.000e-02 9.997e-03 9.998e-03 9.997e-03
9.997e-03 3.000e-02 9.998e-03 9.997e-03
9.998e-03 9.998e-03 3.000e-02 9.998e-03
9.997e-03 9.997e-03 9.998e-03 3.000e-02
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4
1 0.44715 1.000 0.333 0.333 0.333
2 0.44715 0.333 1.000 0.333 0.333
3 0.44717 0.333 0.333 1.000 0.333
4 0.44716 0.333 0.333 0.333 1.000
**********
** 7 **SET ERR 0.5
**********
**********
** 8 **SET PRINT 1
**********
**********
** 9 **HESSE 2000
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=706.56 FROM HESSE STATUS=OK 23 CALLS 89 TOTAL
EDM=1.95881e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 mu_x0 1.80728e-01 1.72984e-01 2.85462e-04 9.04875e-02
2 mu_x1 2.07351e-01 1.72982e-01 2.85769e-04 1.03862e-01
3 mu_x2 -1.59412e-02 1.72987e-01 2.84460e-04 -7.97069e-03
4 mu_x3 1.23430e-01 1.72986e-01 2.84944e-04 6.17540e-02
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 4 ERR DEF=0.5
3.000e-02 9.999e-03 9.999e-03 9.999e-03
9.999e-03 3.000e-02 9.999e-03 9.999e-03
9.999e-03 9.999e-03 3.000e-02 9.999e-03
9.999e-03 9.999e-03 9.999e-03 3.000e-02
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4
1 0.44720 1.000 0.333 0.333 0.333
2 0.44719 0.333 1.000 0.333 0.333
3 0.44720 0.333 0.333 1.000 0.333
4 0.44720 0.333 0.333 0.333 1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_mvg_FOR_OBS_x0:x1:x2:x3 with 0 entries
Metropolis-Hastings progress: ....................................................................................................
[#1] INFO:Eval -- Proposal acceptance rate: 37.1%
[#1] INFO:Eval -- Number of steps in chain: 3710
[#1] INFO:Minization -- createNLL picked up cached constraints from workspace with 0 entries
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --
RooFitResult: minimized FCN value: 706.56, estimated distance to minimum: 1.16082e-08
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
mu_x0 1.8119e-01 +/- 1.73e-01
mu_x1 2.0792e-01 +/- 1.73e-01
mu_x2 -1.6078e-02 +/- 1.73e-01
mu_x3 1.2370e-01 +/- 1.73e-01
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_mvg_mvgData_Profile[mu_x0,mu_x1]) Creating instance of MINUIT
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_mvg_mvgData_Profile[mu_x0,mu_x1]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_mvg_mvgData_Profile[mu_x0,mu_x1]) minimum found at (mu_x0=0.181184, mu_x1=0.207917)
..[#1] INFO:Minization -- LikelihoodInterval - Finding the contour of mu_x0 ( 0 ) and mu_x1 ( 1 )
MCMC interval on p0: [-0.28, 0.6]
MCMC interval on p1: [-0.28, 0.6]
Real time 0:00:01, CP time 1.560
#include "RooGlobalFunc.h"
#include <stdlib.h>
#include "TMatrixDSym.h"
#include "RooArgList.h"
#include "RooRealVar.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "RooAbsReal.h"
#include "RooFitResult.h"
#include "TStopwatch.h"
using namespace std;
using namespace RooFit;
using namespace RooStats;
void MultivariateGaussianTest(Int_t dim = 4, Int_t nPOI = 2)
{
// let's time this challenging example
t.Start();
RooArgList xVec;
RooArgList muVec;
RooArgSet poi;
// make the observable and means
Int_t i, j;
RooRealVar *mu_x;
for (i = 0; i < dim; i++) {
char *name = Form("x%d", i);
x = new RooRealVar(name, name, 0, -3, 3);
char *mu_name = Form("mu_x%d", i);
mu_x = new RooRealVar(mu_name, mu_name, 0, -2, 2);
}
// put them into the list of parameters of interest
for (i = 0; i < nPOI; 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));
ph.SetVariables((RooArgSet &)fit->floatParsFinal());
ph.SetCovMatrix(fit->covarianceMatrix());
ph.SetCacheSize(100);
// 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.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();
}

Definition in file MultivariateGaussianTest.C.

RooMultiVarGaussian.h
RooAbsReal.h
RooStats::MCMCInterval
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: MCMCInterval.h:33
RooStats::ProposalHelper::SetVariables
virtual void SetVariables(RooArgList &vars)
Definition: ProposalHelper.h:54
kGreen
@ kGreen
Definition: Rtypes.h:66
RooStats::ProfileLikelihoodCalculator
The ProfileLikelihoodCalculator is a concrete implementation of CombinedCalculator (the interface cla...
Definition: ProfileLikelihoodCalculator.h:22
Form
char * Form(const char *fmt,...)
TStopwatch.h
RooMultiVarGaussian
Multivariate Gaussian p.d.f.
Definition: RooMultiVarGaussian.h:31
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooStats::ProposalHelper::SetUpdateProposalParameters
virtual void SetUpdateProposalParameters(Bool_t updateParams)
Definition: ProposalHelper.h:51
RooStats::MCMCInterval::LowerLimit
virtual Double_t LowerLimit(RooRealVar &param)
get the lowest value of param that is within the confidence interval
Definition: MCMCInterval.cxx:1006
RooStats::MCMCIntervalPlot::SetLineColor
void SetLineColor(Color_t color)
Definition: MCMCIntervalPlot.h:38
ConfInterval.h
x
Double_t x[n]
Definition: legend1.C:17
TMatrixTSym< Double_t >
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:70
TCanvas.h
RooStats::MCMCIntervalPlot
This class provides simple and straightforward utilities to plot a MCMCInterval object.
Definition: MCMCIntervalPlot.h:28
RooStats::ProposalHelper::SetCacheSize
virtual void SetCacheSize(Int_t size)
Definition: ProposalHelper.h:42
RooFitResult
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
RooStats::ModelConfig::SetPdf
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:81
TMatrixDSym.h
RooStats::LikelihoodIntervalPlot
This class provides simple and straightforward utilities to plot a LikelihoodInterval object.
Definition: LikelihoodIntervalPlot.h:30
TStopwatch::Print
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
RooStats::ProposalHelper
Definition: ProposalHelper.h:29
MarkovChain.h
PyTorch_Generate_CNN_Model.fit
Definition: PyTorch_Generate_CNN_Model.py:32
LikelihoodInterval.h
RooStats::MCMCIntervalPlot::SetLineWidth
void SetLineWidth(Int_t width)
Definition: MCMCIntervalPlot.h:39
MCMCIntervalPlot.h
ModelConfig.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
RooStats::MCMCInterval::UpperLimit
virtual Double_t UpperLimit(RooRealVar &param)
get the highest value of param that is within the confidence interval
Definition: MCMCInterval.cxx:1022
PdfProposal.h
MCMCInterval.h
Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Definition: RooArgSet.cxx:261
RooStats::MCMCCalculator
Bayesian Calculator estimating an interval or a credible region using the Markov-Chain Monte Carlo me...
Definition: MCMCCalculator.h:31
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:437
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
RooRealVar.h
TH2F.h
RooFitResult.h
RooStats::MCMCInterval::GetAxes
virtual RooArgList * GetAxes()
return a list of RooRealVars representing the axes you own the returned RooArgList
Definition: MCMCInterval.h:102
RooGlobalFunc.h
RooStats::ProposalHelper::SetCovMatrix
virtual void SetCovMatrix(const TMatrixDSym &covMatrix)
Definition: ProposalHelper.h:69
MCMCCalculator.h
MetropolisHastings.h
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
Double_t
double Double_t
Definition: RtypesCore.h:59
RooStats::MCMCIntervalPlot::DrawChainScatter
void DrawChainScatter(RooRealVar &xVar, RooRealVar &yVar)
Definition: MCMCIntervalPlot.cxx:590
TCanvas
The Canvas class.
Definition: TCanvas.h:23
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
ProposalFunction.h
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:28
ProfileLikelihoodCalculator.h
name
char name[80]
Definition: TGX11.cxx:110
RooStats::ProposalFunction
ProposalFunction is an interface for all proposal functions that would be used with a Markov Chain Mo...
Definition: ProposalFunction.h:42
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
RooStats::MCMCIntervalPlot::Draw
void Draw(const Option_t *options=NULL)
Definition: MCMCIntervalPlot.cxx:147
RooStats::ModelConfig::SetParametersOfInterest
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
Definition: ModelConfig.h:100
LikelihoodIntervalPlot.h
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
RooFit::Save
RooCmdArg Save(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:190
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
RooArgList.h
RooStats::ProposalHelper::GetProposalFunction
virtual ProposalFunction * GetProposalFunction()
Definition: ProposalHelper.cxx:77
RooStats::LikelihoodInterval
LikelihoodInterval is a concrete implementation of the RooStats::ConfInterval interface.
Definition: LikelihoodInterval.h:34
TCanvas::Update
void Update() override
Definition: TCanvas.cxx:2502
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:182
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
ProposalHelper.h
int
c1
return c1
Definition: legend1.C:41