{
"cells": [
{
"cell_type": "markdown",
"id": "9ab887cb",
"metadata": {},
"source": [
"# MultivariateGaussianTest\n",
"Comparison of MCMC and PLC in a multi-variate gaussian problem\n",
"\n",
"This tutorial produces an N-dimensional multivariate Gaussian\n",
"with a non-trivial covariance matrix. By default N=4 (called \"dim\").\n",
"\n",
"A subset of these are considered parameters of interest.\n",
"This problem is tractable analytically.\n",
"\n",
"We use this mainly as a test of Markov Chain Monte Carlo\n",
"and we compare the result to the profile likelihood ratio.\n",
"\n",
"We use the proposal helper to create a customized\n",
"proposal function for this problem.\n",
"\n",
"For N=4 and 2 parameters of interest it takes about 10-20 seconds\n",
"and the acceptance rate is 37%\n",
"\n",
"\n",
"\n",
"\n",
"**Author:** Artem Busorgin, Kevin Belasco and Kyle Cranmer (C++ version) \n",
"This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, May 19, 2026 at 08:35 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "34d802d8",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:44.001077Z",
"iopub.status.busy": "2026-05-19T20:35:44.000957Z",
"iopub.status.idle": "2026-05-19T20:35:44.969844Z",
"shell.execute_reply": "2026-05-19T20:35:44.969230Z"
}
},
"outputs": [],
"source": [
"import ROOT\n",
"\n",
"dim = 4\n",
"nPOI = 2"
]
},
{
"cell_type": "markdown",
"id": "a4e62351",
"metadata": {},
"source": [
"let's time this challenging example"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "69b57712",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:44.971495Z",
"iopub.status.busy": "2026-05-19T20:35:44.971363Z",
"iopub.status.idle": "2026-05-19T20:35:45.096276Z",
"shell.execute_reply": "2026-05-19T20:35:45.095490Z"
}
},
"outputs": [],
"source": [
"t = ROOT.TStopwatch()\n",
"t.Start()\n",
"\n",
"xVec = []\n",
"muVec = []\n",
"poi = set()"
]
},
{
"cell_type": "markdown",
"id": "88c96d0b",
"metadata": {},
"source": [
"make the observable and means"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "58f40eb2",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.100750Z",
"iopub.status.busy": "2026-05-19T20:35:45.100598Z",
"iopub.status.idle": "2026-05-19T20:35:45.253056Z",
"shell.execute_reply": "2026-05-19T20:35:45.252446Z"
}
},
"outputs": [],
"source": [
"for i in range(dim):\n",
" name = \"x{}\".format(i)\n",
" x = ROOT.RooRealVar(name, name, 0, -3, 3)\n",
" xVec.append(x)\n",
"\n",
" mu_name = \"mu_x{}\".format(i)\n",
" mu_x = ROOT.RooRealVar(mu_name, mu_name, 0, -2, 2)\n",
" muVec.append(mu_x)"
]
},
{
"cell_type": "markdown",
"id": "d73603c7",
"metadata": {},
"source": [
"put them into the list of parameters of interest"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "abc34012",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.256444Z",
"iopub.status.busy": "2026-05-19T20:35:45.256308Z",
"iopub.status.idle": "2026-05-19T20:35:45.364442Z",
"shell.execute_reply": "2026-05-19T20:35:45.363913Z"
}
},
"outputs": [],
"source": [
"for i in range(nPOI):\n",
" poi.add(muVec[i])"
]
},
{
"cell_type": "markdown",
"id": "480c03b7",
"metadata": {},
"source": [
"make a covariance matrix that is all 1's"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "f113e72a",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.366266Z",
"iopub.status.busy": "2026-05-19T20:35:45.366147Z",
"iopub.status.idle": "2026-05-19T20:35:45.516961Z",
"shell.execute_reply": "2026-05-19T20:35:45.516466Z"
}
},
"outputs": [],
"source": [
"cov = ROOT.TMatrixDSym(dim)\n",
"for i in range(dim):\n",
" for j in range(dim):\n",
" if i == j:\n",
" cov[i, j] = 3.0\n",
" else:\n",
" cov[i, j] = 1.0"
]
},
{
"cell_type": "markdown",
"id": "cd3c2a17",
"metadata": {},
"source": [
"now make the multivariate Gaussian"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "256656c7",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.519338Z",
"iopub.status.busy": "2026-05-19T20:35:45.519220Z",
"iopub.status.idle": "2026-05-19T20:35:45.702002Z",
"shell.execute_reply": "2026-05-19T20:35:45.701468Z"
}
},
"outputs": [],
"source": [
"mvg = ROOT.RooMultiVarGaussian(\"mvg\", \"mvg\", xVec, muVec, cov)"
]
},
{
"cell_type": "markdown",
"id": "27cab15f",
"metadata": {},
"source": [
"--------------------\n",
"make a toy dataset"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "d8aed9da",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.703767Z",
"iopub.status.busy": "2026-05-19T20:35:45.703634Z",
"iopub.status.idle": "2026-05-19T20:35:45.836350Z",
"shell.execute_reply": "2026-05-19T20:35:45.835394Z"
}
},
"outputs": [],
"source": [
"data = mvg.generate(xVec, 100)"
]
},
{
"cell_type": "markdown",
"id": "10c844e6",
"metadata": {},
"source": [
"now create the model config for this problem"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "193e3497",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:45.837776Z",
"iopub.status.busy": "2026-05-19T20:35:45.837643Z",
"iopub.status.idle": "2026-05-19T20:35:46.024223Z",
"shell.execute_reply": "2026-05-19T20:35:46.023261Z"
}
},
"outputs": [],
"source": [
"w = ROOT.RooWorkspace(\"MVG\")\n",
"modelConfig = ROOT.RooStats.ModelConfig(w)\n",
"modelConfig.SetPdf(mvg)\n",
"modelConfig.SetParametersOfInterest(poi)"
]
},
{
"cell_type": "markdown",
"id": "118da99d",
"metadata": {},
"source": [
"-------------------------------------------------------\n",
"Setup calculators"
]
},
{
"cell_type": "markdown",
"id": "48be22ba",
"metadata": {},
"source": [
"MCMC\n",
"we want to setup an efficient proposal function\n",
"using the covariance matrix from a fit to the data"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b635c2f8",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:46.025646Z",
"iopub.status.busy": "2026-05-19T20:35:46.025508Z",
"iopub.status.idle": "2026-05-19T20:35:46.278234Z",
"shell.execute_reply": "2026-05-19T20:35:46.273664Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(mvg_over_mvg_Int[x0,x1,x2,x3]) fixing normalization set for coefficient determination to observables in data\n",
"[#1] INFO:Fitting -- using generic CPU library compiled with no vectorizations\n",
"[#1] INFO:Fitting -- Creation of NLL object took 6.7141 ms\n",
"[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_mvg_over_mvg_Int[x0,x1,x2,x3]_mvgData) Summation contains a RooNLLVar, using its error level\n",
"[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"Minuit2Minimizer: Minimize with max-calls 2000 convergence for edm < 1 strategy 1\n",
"Minuit2Minimizer : Valid minimum - status = 0\n",
"FVAL = 706.560063684865781\n",
"Edm = 9.79361268420962158e-06\n",
"Nfcn = 68\n",
"mu_x0\t = 0.180728\t +/- 0.17298\t(limited)\n",
"mu_x1\t = 0.207351\t +/- 0.172978\t(limited)\n",
"mu_x2\t = -0.0159412\t +/- 0.172984\t(limited)\n",
"mu_x3\t = 0.12343\t +/- 0.172982\t(limited)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Info in : MnSeedGenerator Computing seed using NumericalGradient calculator\n",
"Info in : MnSeedGenerator Evaluated function and gradient in 295.671 μs\n",
"Info in : MnSeedGenerator Initial state: FCN = 707.8223914 Edm = 1.103554631 NCalls = 17\n",
"Info in : MnSeedGenerator Initial state \n",
" Minimum value : 707.8223914\n",
" Edm : 1.103554631\n",
" Internal parameters:\t[ 0 0 0 0]\t\n",
" Internal gradient :\t[ -9.839418409 -12.51273889 9.885731533 -4.091554467]\t\n",
" Internal covariance matrix:\n",
"[[ 0.012000008 0 0 0]\n",
" [ 0 0.012000008 0 0]\n",
" [ 0 0 0.012000008 0]\n",
" [ 0 0 0 0.012000008]]]\n",
"Info in : VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 2000\n",
"Info in : VariableMetricBuilder 0 - FCN = 707.8223914 Edm = 1.103554631 NCalls = 17\n",
"Info in : VariableMetricBuilder 1 - FCN = 706.7749736 Edm = 0.1093739779 NCalls = 26\n",
"Info in : VariableMetricBuilder 2 - FCN = 706.5648985 Edm = 0.004643182323 NCalls = 36\n",
"Info in : VariableMetricBuilder 3 - FCN = 706.5600637 Edm = 9.358845203e-06 NCalls = 45\n",
"Info in : MnHesse Done after 273.272 μs\n",
"Info in : VariableMetricBuilder After Hessian\n",
"Info in : VariableMetricBuilder 4 - FCN = 706.5600637 Edm = 9.793612684e-06 NCalls = 68\n",
"Info in : VariableMetricBuilder Stop iterating after 686.884 μs\n",
"Info in : Minuit2Minimizer::Hesse Using max-calls 2000\n",
"Info in : MnHesse Done after 237.121 μs\n",
"Info in : Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate\n"
]
}
],
"source": [
"fit = mvg.fitTo(data, Save=True)\n",
"ph = ROOT.RooStats.ProposalHelper()\n",
"ph.SetVariables(fit.floatParsFinal())\n",
"ph.SetCovMatrix(fit.covarianceMatrix())\n",
"ph.SetUpdateProposalParameters(True)\n",
"ph.SetCacheSize(100)\n",
"pdfProp = ph.GetProposalFunction()"
]
},
{
"cell_type": "markdown",
"id": "1ce5e86c",
"metadata": {},
"source": [
"now create the calculator"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "eb52a5cf",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:46.280067Z",
"iopub.status.busy": "2026-05-19T20:35:46.279924Z",
"iopub.status.idle": "2026-05-19T20:35:46.885136Z",
"shell.execute_reply": "2026-05-19T20:35:46.884460Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(mvg_over_mvg_Int[x0,x1,x2,x3]) fixing normalization set for coefficient determination to observables in data\n",
"[#1] INFO:Fitting -- Creation of NLL object took 176.521 μs\n",
"Metropolis-Hastings progress: ....................................................................................................\n",
"[#1] INFO:Eval -- Proposal acceptance rate: 37.1%\n",
"[#1] INFO:Eval -- Number of steps in chain: 3710\n"
]
}
],
"source": [
"mc = ROOT.RooStats.MCMCCalculator(data, modelConfig)\n",
"mc.SetConfidenceLevel(0.95)\n",
"mc.SetNumBurnInSteps(100)\n",
"mc.SetNumIters(10000)\n",
"mc.SetNumBins(50)\n",
"mc.SetProposalFunction(pdfProp)\n",
"\n",
"mcInt = mc.GetInterval()\n",
"poiList = mcInt.GetAxes()"
]
},
{
"cell_type": "markdown",
"id": "04a998ae",
"metadata": {},
"source": [
"now setup the profile likelihood calculator"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "ba53bce3",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:46.897220Z",
"iopub.status.busy": "2026-05-19T20:35:46.897082Z",
"iopub.status.idle": "2026-05-19T20:35:47.023907Z",
"shell.execute_reply": "2026-05-19T20:35:47.023256Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:InputArguments -- The deprecated RooFit::CloneData(1) option passed to createNLL() is ignored.\n",
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(mvg_over_mvg_Int[x0,x1,x2,x3]) fixing normalization set for coefficient determination to observables in data\n",
"[#1] INFO:Fitting -- Creation of NLL object took 177.309 μs\n",
"[#0] PROGRESS:Minimization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE \n",
"[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_mvg_over_mvg_Int[x0,x1,x2,x3]_mvgData) Summation contains a RooNLLVar, using its error level\n",
"[#0] PROGRESS:Minimization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit2 / Migrad with strategy 1\n",
"[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"[#1] INFO:Minimization -- \n",
" RooFitResult: minimized FCN value: 706.56, estimated distance to minimum: 2.41067e-06\n",
" covariance matrix quality: Full, accurate covariance matrix\n",
" Status : MINIMIZE=0 \n",
"\n",
" Floating Parameter FinalValue +/- Error \n",
" -------------------- --------------------------\n",
" mu_x0 1.8134e-01 +/- 1.73e-01\n",
" mu_x1 2.0772e-01 +/- 1.73e-01\n",
" mu_x2 -1.5923e-02 +/- 1.73e-01\n",
" mu_x3 1.2360e-01 +/- 1.73e-01\n",
"\n"
]
}
],
"source": [
"plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, modelConfig)\n",
"plc.SetConfidenceLevel(0.95)\n",
"plInt = plc.GetInterval()"
]
},
{
"cell_type": "markdown",
"id": "852cee9f",
"metadata": {},
"source": [
"make some plots"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "3175a776",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:47.047232Z",
"iopub.status.busy": "2026-05-19T20:35:47.047082Z",
"iopub.status.idle": "2026-05-19T20:35:47.422545Z",
"shell.execute_reply": "2026-05-19T20:35:47.422207Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[mu_x0,mu_x1]) Creating instance of MINUIT\n",
"[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_mvg_over_mvg_Int[x0,x1,x2,x3]_mvgData) Summation contains a RooNLLVar, using its error level\n",
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[mu_x0,mu_x1]) determining minimum likelihood for current configurations w.r.t all observable\n",
"[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[mu_x0,mu_x1]) minimum found at (mu_x0=0.181185, mu_x1=0.207918)\n",
".[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
".[#1] INFO:Minimization -- [fitFCN] No discrete parameters, performing continuous minimization only\n",
"[#1] INFO:Minimization -- LikelihoodInterval - Finding the contour of mu_x0 ( 0 ) and mu_x1 ( 1 ) \n",
"MCMC interval on p0: [-0.27999999999999997, 0.6000000000000001]\n",
"MCMC interval on p1: [-0.19999999999999996, 0.6000000000000001]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Real time 0:00:02, CP time 1.060\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Info in : png file MultivariateGaussianTest_scatter.png has been created\n",
"Info in : png file MultivariateGaussianTest_plot.png has been created\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mcPlot = ROOT.RooStats.MCMCIntervalPlot(mcInt)\n",
"\n",
"c1 = ROOT.TCanvas()\n",
"mcPlot.SetLineColor(\"kGreen\")\n",
"mcPlot.SetLineWidth(2)\n",
"mcPlot.Draw()\n",
"\n",
"plPlot = ROOT.RooStats.LikelihoodIntervalPlot(plInt)\n",
"plPlot.Draw(\"same\")\n",
"\n",
"if len(poiList) == 1:\n",
" p = poiList.at(0)\n",
" print(\"MCMC interval: [{}, {}]\".format(mcInt.LowerLimit(p), mcInt.UpperLimit(p)))\n",
"\n",
"if len(poiList) == 2:\n",
" p0 = poiList.at(0)\n",
" p1 = poiList.at(1)\n",
" scatter = ROOT.TCanvas()\n",
" print(\"MCMC interval on p0: [{}, {}]\".format(mcInt.LowerLimit(p0), mcInt.UpperLimit(p0)))\n",
" print(\"MCMC interval on p1: [{}, {}]\".format(mcInt.LowerLimit(p1), mcInt.UpperLimit(p1)))\n",
"\n",
" # MCMC interval on p0: [-0.2, 0.6]\n",
" # MCMC interval on p1: [-0.2, 0.6]\n",
"\n",
" mcPlot.DrawChainScatter(p0, p1)\n",
" scatter.Update()\n",
" scatter.SaveAs(\"MultivariateGaussianTest_scatter.png\")\n",
"\n",
"t.Stop()\n",
"t.Print()\n",
"\n",
"c1.SaveAs(\"MultivariateGaussianTest_plot.png\")"
]
},
{
"cell_type": "markdown",
"id": "f592205a",
"metadata": {},
"source": [
"TODO: The MCMCCalculator has to be destructed first. Otherwise, we can get\n",
"segmentation faults depending on the destruction order, which is random in\n",
"Python. Probably the issue is that some object has a non-owning pointer to\n",
"another object, which it uses in its destructor. This should be fixed either\n",
"in the design of RooStats in C++, or with phythonizations."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "23e1e348",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:47.440264Z",
"iopub.status.busy": "2026-05-19T20:35:47.440121Z",
"iopub.status.idle": "2026-05-19T20:35:47.550205Z",
"shell.execute_reply": "2026-05-19T20:35:47.549626Z"
}
},
"outputs": [],
"source": [
"del mc"
]
},
{
"cell_type": "markdown",
"id": "26337059",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "087a4384",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2026-05-19T20:35:47.552075Z",
"iopub.status.busy": "2026-05-19T20:35:47.551924Z",
"iopub.status.idle": "2026-05-19T20:35:47.673423Z",
"shell.execute_reply": "2026-05-19T20:35:47.672731Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from ROOT import gROOT \n",
"gROOT.GetListOfCanvases().Draw()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}