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%
Processing /mnt/build/workspace/root-makedoc-v610/rootspi/rdoc/src/v6-10-00-patches/tutorials/roostats/MultivariateGaussianTest.C...
#include <stdlib.h>
void MultivariateGaussianTest(
Int_t dim = 4,
Int_t nPOI = 2)
{
for (i = 0; i < dim; i++) {
char* name =
Form(
"x%d", i);
char* mu_name =
Form(
"mu_x%d",i);
}
for (i = 0; i < nPOI; i++) {
}
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;
}
}
mc.SetConfidenceLevel(0.95);
mc.SetNumBurnInSteps(100);
mc.SetNumIters(10000);
mc.SetNumBins(50);
mc.SetProposalFunction(*pdfProp);
plc.SetConfidenceLevel(0.95);
plPlot.Draw("same");
cout << "MCMC interval: [" << ll << ", " << ul << "]" << endl;
}
cout << "MCMC interval on p0: [" << ll << ", " << ul << "]" << endl;
cout << "MCMC interval on p1: [" << ll << ", " << ul << "]" << endl;
}
}
- Authors
- Kevin Belasco and Kyle Cranmer
Definition in file MultivariateGaussianTest.C.