{
"cells": [
{
"cell_type": "markdown",
"id": "7ef51b55",
"metadata": {},
"source": [
"# rs302_JeffreysPriorDemo\n",
"tutorial demonstrating and validates the RooJeffreysPrior class\n",
"\n",
"Jeffreys's prior is an 'objective prior' based on formal rules.\n",
"It is calculated from the Fisher information matrix.\n",
"\n",
"Read more:\n",
"http://en.wikipedia.org/wiki/Jeffreys_prior\n",
"\n",
"The analytic form is not known for most PDFs, but it is for\n",
"simple cases like the Poisson mean, Gaussian mean, Gaussian sigma.\n",
"\n",
"This class uses numerical tricks to calculate the Fisher Information Matrix\n",
"efficiently. In particular, it takes advantage of a property of the\n",
"'Asimov data' as described in\n",
"Asymptotic formulae for likelihood-based tests of new physics\n",
"Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells\n",
"http://arxiv.org/abs/arXiv:1007.1727\n",
"\n",
"This Demo has four parts:\n",
" 1. TestJeffreysPriorDemo -- validates Poisson mean case 1/sqrt(mu)\n",
" 2. TestJeffreysGaussMean -- validates Gaussian mean case\n",
" 3. TestJeffreysGaussSigma -- validates Gaussian sigma case 1/sigma\n",
" 4. TestJeffreysGaussMeanAndSigma -- demonstrates 2-d example\n",
"\n",
"\n",
"\n",
"\n",
"**Author:** Kyle Cranmer \n",
"This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:18 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "4aabdf04",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:29.434130Z",
"iopub.status.busy": "2024-03-19T19:18:29.433750Z",
"iopub.status.idle": "2024-03-19T19:18:29.470275Z",
"shell.execute_reply": "2024-03-19T19:18:29.469220Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"#include \"RooJeffreysPrior.h\"\n",
"\n",
"#include \"RooWorkspace.h\"\n",
"#include \"RooDataHist.h\"\n",
"#include \"RooGenericPdf.h\"\n",
"#include \"RooPlot.h\"\n",
"#include \"RooFitResult.h\"\n",
"#include \"RooRealVar.h\"\n",
"#include \"RooAbsPdf.h\"\n",
"#include \"RooNumIntConfig.h\"\n",
"\n",
"#include \"TH1F.h\"\n",
"#include \"TCanvas.h\"\n",
"#include \"TLegend.h\"\n",
"#include \"TMatrixDSym.h\"\n",
"\n",
"using namespace RooFit;"
]
},
{
"cell_type": "markdown",
"id": "c9e6a327",
"metadata": {},
"source": [
" _________________________________________________\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "99adf69f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:29.476793Z",
"iopub.status.busy": "2024-03-19T19:18:29.476384Z",
"iopub.status.idle": "2024-03-19T19:18:29.723529Z",
"shell.execute_reply": "2024-03-19T19:18:29.722191Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"void TestJeffreysGaussMean()\n",
"{\n",
" RooWorkspace w(\"w\");\n",
" w.factory(\"Gaussian::g(x[0,-20,20],mu[0,-5.,5],sigma[1,0,10])\");\n",
" w.factory(\"n[10,.1,200]\");\n",
" w.factory(\"ExtendPdf::p(g,n)\");\n",
" w.var(\"sigma\")->setConstant();\n",
" w.var(\"n\")->setConstant();\n",
"\n",
" std::unique_ptr asimov{w.pdf(\"p\")->generateBinned(*w.var(\"x\"), ExpectedData())};\n",
"\n",
" std::unique_ptr res{w.pdf(\"p\")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};\n",
"\n",
" asimov->Print();\n",
" res->Print();\n",
" TMatrixDSym cov = res->covarianceMatrix();\n",
" cout << \"variance = \" << (cov.Determinant()) << endl;\n",
" cout << \"stdev = \" << sqrt(cov.Determinant()) << endl;\n",
" cov.Invert();\n",
" cout << \"jeffreys = \" << sqrt(cov.Determinant()) << endl;\n",
"\n",
" w.defineSet(\"poi\", \"mu\");\n",
" w.defineSet(\"obs\", \"x\");\n",
"\n",
" RooJeffreysPrior pi(\"jeffreys\", \"jeffreys\", *w.pdf(\"p\"), *w.set(\"poi\"), *w.set(\"obs\"));\n",
"\n",
" const RooArgSet *temp = w.set(\"poi\");\n",
" pi.getParameters(*temp)->Print();\n",
"\n",
" // return;\n",
" RooGenericPdf *test = new RooGenericPdf(\"constant\", \"Expected = constant\", \"1\", *w.set(\"poi\"));\n",
"\n",
" TCanvas *c1 = new TCanvas;\n",
" RooPlot *plot = w.var(\"mu\")->frame();\n",
" pi.plotOn(plot);\n",
" test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));\n",
" plot->Draw();\n",
"\n",
" auto legend = plot->BuildLegend();\n",
" legend->DrawClone();\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "bea1862a",
"metadata": {},
"source": [
" _________________________________________________\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3a61b957",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:29.728518Z",
"iopub.status.busy": "2024-03-19T19:18:29.728142Z",
"iopub.status.idle": "2024-03-19T19:18:29.828450Z",
"shell.execute_reply": "2024-03-19T19:18:29.821154Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"void TestJeffreysGaussSigma()\n",
"{\n",
" // this one is VERY sensitive\n",
" // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved\n",
" // and you get really bizarre shapes\n",
" // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized\n",
" // and the PDF falls off too fast at high sigma\n",
" RooWorkspace w(\"w\");\n",
" w.factory(\"Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])\");\n",
" w.factory(\"n[100,.1,2000]\");\n",
" w.factory(\"ExtendPdf::p(g,n)\");\n",
" // w.var(\"sigma\")->setConstant();\n",
" w.var(\"mu\")->setConstant();\n",
" w.var(\"n\")->setConstant();\n",
" w.var(\"x\")->setBins(301);\n",
"\n",
" std::unique_ptr asimov{w.pdf(\"p\")->generateBinned(*w.var(\"x\"), ExpectedData())};\n",
"\n",
" std::unique_ptr res{w.pdf(\"p\")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};\n",
"\n",
" asimov->Print();\n",
" res->Print();\n",
" TMatrixDSym cov = res->covarianceMatrix();\n",
" cout << \"variance = \" << (cov.Determinant()) << endl;\n",
" cout << \"stdev = \" << sqrt(cov.Determinant()) << endl;\n",
" cov.Invert();\n",
" cout << \"jeffreys = \" << sqrt(cov.Determinant()) << endl;\n",
"\n",
" w.defineSet(\"poi\", \"sigma\");\n",
" w.defineSet(\"obs\", \"x\");\n",
"\n",
" RooJeffreysPrior pi(\"jeffreys\", \"jeffreys\", *w.pdf(\"p\"), *w.set(\"poi\"), *w.set(\"obs\"));\n",
"\n",
" const RooArgSet *temp = w.set(\"poi\");\n",
" pi.getParameters(*temp)->Print();\n",
"\n",
" RooGenericPdf *test = new RooGenericPdf(\"test\", \"Expected = #sqrt{2}/#sigma\", \"sqrt(2.)/sigma\", *w.set(\"poi\"));\n",
"\n",
" TCanvas *c1 = new TCanvas;\n",
" RooPlot *plot = w.var(\"sigma\")->frame();\n",
" pi.plotOn(plot);\n",
" test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));\n",
" plot->Draw();\n",
"\n",
" auto legend = plot->BuildLegend();\n",
" legend->DrawClone();\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "696674cd",
"metadata": {},
"source": [
" _________________________________________________\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3380b438",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:29.846585Z",
"iopub.status.busy": "2024-03-19T19:18:29.846195Z",
"iopub.status.idle": "2024-03-19T19:18:29.928828Z",
"shell.execute_reply": "2024-03-19T19:18:29.927395Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"void TestJeffreysGaussMeanAndSigma()\n",
"{\n",
" // this one is VERY sensitive\n",
" // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved\n",
" // and you get really bizarre shapes\n",
" // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized\n",
" // and the PDF falls off too fast at high sigma\n",
" RooWorkspace w(\"w\");\n",
" w.factory(\"Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1.,5.])\");\n",
" w.factory(\"n[100,.1,2000]\");\n",
" w.factory(\"ExtendPdf::p(g,n)\");\n",
"\n",
" w.var(\"n\")->setConstant();\n",
" w.var(\"x\")->setBins(301);\n",
"\n",
" std::unique_ptr asimov{w.pdf(\"p\")->generateBinned(*w.var(\"x\"), ExpectedData())};\n",
"\n",
" std::unique_ptr res{w.pdf(\"p\")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};\n",
"\n",
" asimov->Print();\n",
" res->Print();\n",
" TMatrixDSym cov = res->covarianceMatrix();\n",
" cout << \"variance = \" << (cov.Determinant()) << endl;\n",
" cout << \"stdev = \" << sqrt(cov.Determinant()) << endl;\n",
" cov.Invert();\n",
" cout << \"jeffreys = \" << sqrt(cov.Determinant()) << endl;\n",
"\n",
" w.defineSet(\"poi\", \"mu,sigma\");\n",
" w.defineSet(\"obs\", \"x\");\n",
"\n",
" RooJeffreysPrior pi(\"jeffreys\", \"jeffreys\", *w.pdf(\"p\"), *w.set(\"poi\"), *w.set(\"obs\"));\n",
"\n",
" const RooArgSet *temp = w.set(\"poi\");\n",
" pi.getParameters(*temp)->Print();\n",
" // return;\n",
"\n",
" TCanvas *c1 = new TCanvas;\n",
" TH1 *Jeff2d = pi.createHistogram(\"2dJeffreys\", *w.var(\"mu\"), Binning(10, -5., 5), YVar(*w.var(\"sigma\"), Binning(10, 1., 5.)));\n",
" Jeff2d->Draw(\"surf\");\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "2aa769a3",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:29.938327Z",
"iopub.status.busy": "2024-03-19T19:18:29.937914Z",
"iopub.status.idle": "2024-03-19T19:18:31.805883Z",
"shell.execute_reply": "2024-03-19T19:18:31.804779Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.\n",
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(p) fixing normalization set for coefficient determination to observables in data\n",
"[#1] INFO:Fitting -- using CPU computation library compiled with -mavx2\n",
"[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_p_genData) Summation contains a RooNLLVar, using its error level\n",
"[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization\n",
"Minuit2Minimizer: Minimize with max-calls 500 convergence for edm < 1 strategy 1\n",
"Minuit2Minimizer : Valid minimum - status = 0\n",
"FVAL = -360.517018598809159\n",
"Edm = 4.13877261906962124e-17\n",
"Nfcn = 22\n",
"mu\t = 100\t +/- 9.98317\t(limited)\n",
"[#1] INFO:Fitting -- RooAbsPdf::fitTo(p) Calculating sum-of-weights-squared correction matrix for covariance matrix\n",
"[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization\n",
"RooDataHist::genData[x] = 100 bins (100 weights)\n",
"\n",
" RooFitResult: minimized FCN value: -360.517, estimated distance to minimum: 3.66543e-18\n",
" covariance matrix quality: Full, accurate covariance matrix\n",
" Status : MINIMIZE=0 HESSE=0 HESSE=0 \n",
"\n",
" Floating Parameter FinalValue +/- Error \n",
" -------------------- --------------------------\n",
" mu 1.0000e+02 +/- 1.00e+01\n",
"\n",
"variance = 100.016\n",
"stdev = 10.0008\n",
"jeffreys = 0.0999921\n",
"[#1] INFO:NumericIntegration -- RooRealIntegral::init(jeffreys_Int[mu]) using numeric integrator RooAdaptiveGaussKronrodIntegrator1D to calculate Int(mu)\n",
"[#1] INFO:NumericIntegration -- RooRealIntegral::init(Expected_Int[mu]) using numeric integrator RooIntegrator1D to calculate Int(mu)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Info in : MnSeedGenerator Computing seed using NumericalGradient calculator\n",
"Info in : MnSeedGenerator Initial state: FCN = -360.5170186 Edm = 1.621500752e-11 NCalls = 5\n",
"Info in : MnSeedGenerator Initial state \n",
" Minimum value : -360.5170186\n",
" Edm : 1.621500752e-11\n",
" Internal parameters:\t[ -0.005025146777]\t\n",
" Internal gradient :\t[ -5.666191281e-05]\t\n",
" Internal covariance matrix:\n",
"[[ 0.020202015]]]\n",
"Info in : VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 500\n",
"Info in : VariableMetricBuilder 0 - FCN = -360.5170186 Edm = 1.621500752e-11 NCalls = 5\n",
"Warning in : VariableMetricBuilder No improvement in line search\n",
"Info in : VariableMetricBuilder 1 - FCN = -360.5170186 Edm = 1.621500752e-11 NCalls = 15\n",
"Info in : VariableMetricBuilder After Hessian\n",
"Info in : VariableMetricBuilder 2 - FCN = -360.5170186 Edm = 4.138772619e-17 NCalls = 22\n",
"Info in : Minuit2Minimizer::Hesse Using max-calls 500\n",
"Info in : Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate\n",
"Info in : Minuit2Minimizer::Hesse Using max-calls 500\n",
"Info in : Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate\n"
]
}
],
"source": [
"RooWorkspace w(\"w\");\n",
"w.factory(\"Uniform::u(x[0,1])\");\n",
"w.factory(\"mu[100,1,200]\");\n",
"w.factory(\"ExtendPdf::p(u,mu)\");\n",
"\n",
"std::unique_ptr asimov{w.pdf(\"p\")->generateBinned(*w.var(\"x\"), ExpectedData())};\n",
"\n",
"std::unique_ptr res{w.pdf(\"p\")->fitTo(*asimov, Save(), SumW2Error(kTRUE))};\n",
"\n",
"asimov->Print();\n",
"res->Print();\n",
"TMatrixDSym cov = res->covarianceMatrix();\n",
"cout << \"variance = \" << (cov.Determinant()) << endl;\n",
"cout << \"stdev = \" << sqrt(cov.Determinant()) << endl;\n",
"cov.Invert();\n",
"cout << \"jeffreys = \" << sqrt(cov.Determinant()) << endl;\n",
"\n",
"w.defineSet(\"poi\", \"mu\");\n",
"w.defineSet(\"obs\", \"x\");\n",
"\n",
"RooJeffreysPrior pi(\"jeffreys\", \"jeffreys\", *w.pdf(\"p\"), *w.set(\"poi\"), *w.set(\"obs\"));\n",
"\n",
"RooGenericPdf *test = new RooGenericPdf(\"Expected\", \"Expected = 1/#sqrt{#mu}\", \"1./sqrt(mu)\",\n",
" *w.set(\"poi\"));\n",
"\n",
"TCanvas *c1 = new TCanvas;\n",
"RooPlot *plot = w.var(\"mu\")->frame();\n",
"pi.plotOn(plot);\n",
"test->plotOn(plot, LineColor(kRed), LineStyle(kDashDotted));\n",
"plot->Draw();\n",
"\n",
"auto legend = plot->BuildLegend();\n",
"legend->DrawClone();"
]
},
{
"cell_type": "markdown",
"id": "4456ee0c",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "2af466f3",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:31.817925Z",
"iopub.status.busy": "2024-03-19T19:18:31.817535Z",
"iopub.status.idle": "2024-03-19T19:18:32.218735Z",
"shell.execute_reply": "2024-03-19T19:18:32.217122Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%jsroot on\n",
"gROOT->GetListOfCanvases()->Draw()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "ROOT C++",
"language": "c++",
"name": "root"
},
"language_info": {
"codemirror_mode": "text/x-c++src",
"file_extension": ".C",
"mimetype": " text/x-c++src",
"name": "c++"
}
},
"nbformat": 4,
"nbformat_minor": 5
}