{ "cells": [ { "cell_type": "markdown", "id": "9c86ce6f", "metadata": {}, "source": [ "# IntervalExamples\n", "Example showing confidence intervals with four techniques.\n", "\n", "An example that shows confidence intervals with four techniques.\n", "The model is a Normal Gaussian G(x|mu,sigma) with 100 samples of x.\n", "The answer is known analytically, so this is a good example to validate\n", "the RooStats tools.\n", "\n", " - expected interval is [-0.162917, 0.229075]\n", " - plc interval is [-0.162917, 0.229075]\n", " - fc interval is [-0.17 , 0.23] // stepsize is 0.01\n", " - bc interval is [-0.162918, 0.229076]\n", " - mcmc interval is [-0.166999, 0.230224]\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, May 19, 2026 at 08:35 PM." ] }, { "cell_type": "code", "execution_count": null, "id": "add3bd93", "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cpp -d\n", "#include \"RooStats/ConfInterval.h\"\n", "#include \"RooStats/PointSetInterval.h\"\n", "#include \"RooStats/ConfidenceBelt.h\"\n", "#include \"RooStats/FeldmanCousins.h\"\n", "#include \"RooStats/ProfileLikelihoodCalculator.h\"\n", "#include \"RooStats/MCMCCalculator.h\"\n", "#include \"RooStats/BayesianCalculator.h\"\n", "#include \"RooStats/MCMCIntervalPlot.h\"\n", "#include \"RooStats/LikelihoodIntervalPlot.h\"\n", "\n", "#include \"RooStats/ToyMCSampler.h\"\n", "\n", "#include \"RooRandom.h\"\n", "#include \"RooDataSet.h\"\n", "#include \"RooRealVar.h\"\n", "#include \"RooConstVar.h\"\n", "#include \"RooAddition.h\"\n", "#include \"RooDataHist.h\"\n", "#include \"RooPoisson.h\"\n", "#include \"RooPlot.h\"\n", "\n", "#include \"TCanvas.h\"\n", "#include \"TTree.h\"\n", "#include \"TStyle.h\"\n", "#include \"TMath.h\"\n", "#include \"Math/DistFunc.h\"\n", "#include \"TH1F.h\"\n", "#include \"TMarker.h\"\n", "#include \"TStopwatch.h\"\n", "\n", "#include \n", "\n", "using namespace RooFit;\n", "using namespace RooStats;" ] }, { "cell_type": "markdown", "id": "9ac9c59d", "metadata": {}, "source": [ "Time this macro" ] }, { "cell_type": "code", "execution_count": null, "id": "56d736fa", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TStopwatch t;\n", "t.Start();" ] }, { "cell_type": "markdown", "id": "fb4d1eb1", "metadata": {}, "source": [ "set RooFit random seed for reproducible results" ] }, { "cell_type": "code", "execution_count": null, "id": "2ca3efb4", "metadata": { "collapsed": false }, "outputs": [], "source": [ "RooRandom::randomGenerator()->SetSeed(3001);" ] }, { "cell_type": "markdown", "id": "95be378e", "metadata": {}, "source": [ "make a simple model via the workspace factory" ] }, { "cell_type": "code", "execution_count": null, "id": "bf401c35", "metadata": { "collapsed": false }, "outputs": [], "source": [ "RooWorkspace *wspace = new RooWorkspace();\n", "wspace->factory(\"Gaussian::normal(x[-10,10],mu[-1,1],sigma[1])\");\n", "wspace->defineSet(\"poi\", \"mu\");\n", "wspace->defineSet(\"obs\", \"x\");" ] }, { "cell_type": "markdown", "id": "be44e46d", "metadata": {}, "source": [ "specify components of model for statistical tools" ] }, { "cell_type": "code", "execution_count": null, "id": "6bc4d514", "metadata": { "collapsed": false }, "outputs": [], "source": [ "ModelConfig *modelConfig = new ModelConfig(\"Example G(x|mu,1)\");\n", "modelConfig->SetWorkspace(*wspace);\n", "modelConfig->SetPdf(*wspace->pdf(\"normal\"));\n", "modelConfig->SetParametersOfInterest(*wspace->set(\"poi\"));\n", "modelConfig->SetObservables(*wspace->set(\"obs\"));" ] }, { "cell_type": "markdown", "id": "089b2e1e", "metadata": {}, "source": [ "create a toy dataset" ] }, { "cell_type": "code", "execution_count": null, "id": "bf5db2dc", "metadata": { "collapsed": false }, "outputs": [], "source": [ "std::unique_ptr data{wspace->pdf(\"normal\")->generate(*wspace->set(\"obs\"), 100)};\n", "data->Print();" ] }, { "cell_type": "markdown", "id": "91bd1eed", "metadata": {}, "source": [ "for convenience later on" ] }, { "cell_type": "code", "execution_count": null, "id": "1d1eaad3", "metadata": { "collapsed": false }, "outputs": [], "source": [ "RooRealVar *x = wspace->var(\"x\");\n", "RooRealVar *mu = wspace->var(\"mu\");" ] }, { "cell_type": "markdown", "id": "5bc94e54", "metadata": {}, "source": [ "set confidence level" ] }, { "cell_type": "code", "execution_count": null, "id": "2f5920ae", "metadata": { "collapsed": false }, "outputs": [], "source": [ "double confidenceLevel = 0.95;" ] }, { "cell_type": "markdown", "id": "fd038260", "metadata": {}, "source": [ "example use profile likelihood calculator" ] }, { "cell_type": "code", "execution_count": null, "id": "571ec5e5", "metadata": { "collapsed": false }, "outputs": [], "source": [ "ProfileLikelihoodCalculator plc(*data, *modelConfig);\n", "plc.SetConfidenceLevel(confidenceLevel);\n", "LikelihoodInterval *plInt = plc.GetInterval();" ] }, { "cell_type": "markdown", "id": "f63266b0", "metadata": {}, "source": [ "example use of Feldman-Cousins" ] }, { "cell_type": "code", "execution_count": null, "id": "0ff5db3a", "metadata": { "collapsed": false }, "outputs": [], "source": [ "FeldmanCousins fc(*data, *modelConfig);\n", "fc.SetConfidenceLevel(confidenceLevel);\n", "fc.SetNBins(100); // number of points to test per parameter\n", "fc.UseAdaptiveSampling(true); // make it go faster" ] }, { "cell_type": "markdown", "id": "48e9842d", "metadata": {}, "source": [ "Here, we consider only ensembles with 100 events\n", "The PDF could be extended and this could be removed" ] }, { "cell_type": "code", "execution_count": null, "id": "c8291830", "metadata": { "collapsed": false }, "outputs": [], "source": [ "fc.FluctuateNumDataEntries(false);\n", "\n", "PointSetInterval *interval = (PointSetInterval *)fc.GetInterval();" ] }, { "cell_type": "markdown", "id": "726f1762", "metadata": {}, "source": [ "example use of BayesianCalculator\n", "now we also need to specify a prior in the ModelConfig" ] }, { "cell_type": "code", "execution_count": null, "id": "a2422656", "metadata": { "collapsed": false }, "outputs": [], "source": [ "wspace->factory(\"Uniform::prior(mu)\");\n", "modelConfig->SetPriorPdf(*wspace->pdf(\"prior\"));" ] }, { "cell_type": "markdown", "id": "ad70f93a", "metadata": {}, "source": [ "example usage of BayesianCalculator" ] }, { "cell_type": "code", "execution_count": null, "id": "167cea0b", "metadata": { "collapsed": false }, "outputs": [], "source": [ "BayesianCalculator bc(*data, *modelConfig);\n", "bc.SetConfidenceLevel(confidenceLevel);\n", "SimpleInterval *bcInt = bc.GetInterval();" ] }, { "cell_type": "markdown", "id": "c4d2830c", "metadata": {}, "source": [ "example use of MCMCInterval" ] }, { "cell_type": "code", "execution_count": null, "id": "08cf0317", "metadata": { "collapsed": false }, "outputs": [], "source": [ "MCMCCalculator mc(*data, *modelConfig);\n", "mc.SetConfidenceLevel(confidenceLevel);" ] }, { "cell_type": "markdown", "id": "c41b6299", "metadata": {}, "source": [ "special options" ] }, { "cell_type": "code", "execution_count": null, "id": "fad8ccdf", "metadata": { "collapsed": false }, "outputs": [], "source": [ "mc.SetNumBins(200); // bins used internally for representing posterior\n", "mc.SetNumBurnInSteps(500); // first N steps to be ignored as burn-in\n", "mc.SetNumIters(100000); // how long to run chain\n", "mc.SetLeftSideTailFraction(0.5); // for central interval\n", "MCMCInterval *mcInt = mc.GetInterval();" ] }, { "cell_type": "markdown", "id": "45666f45", "metadata": {}, "source": [ "for this example we know the expected intervals" ] }, { "cell_type": "code", "execution_count": null, "id": "65c9ee24", "metadata": { "collapsed": false }, "outputs": [], "source": [ "double expectedLL =\n", " data->mean(*x) + ROOT::Math::normal_quantile((1 - confidenceLevel) / 2, 1) / sqrt(data->numEntries());\n", "double expectedUL =\n", " data->mean(*x) + ROOT::Math::normal_quantile_c((1 - confidenceLevel) / 2, 1) / sqrt(data->numEntries());" ] }, { "cell_type": "markdown", "id": "3a9a145a", "metadata": {}, "source": [ "Use the intervals" ] }, { "cell_type": "code", "execution_count": null, "id": "36e3b854", "metadata": { "collapsed": false }, "outputs": [], "source": [ "std::cout << \"expected interval is [\" << expectedLL << \", \" << expectedUL << \"]\" << endl;\n", "\n", "cout << \"plc interval is [\" << plInt->LowerLimit(*mu) << \", \" << plInt->UpperLimit(*mu) << \"]\" << endl;\n", "\n", "std::cout << \"fc interval is [\" << interval->LowerLimit(*mu) << \" , \" << interval->UpperLimit(*mu) << \"]\" << endl;\n", "\n", "cout << \"bc interval is [\" << bcInt->LowerLimit() << \", \" << bcInt->UpperLimit() << \"]\" << endl;\n", "\n", "cout << \"mc interval is [\" << mcInt->LowerLimit(*mu) << \", \" << mcInt->UpperLimit(*mu) << \"]\" << endl;\n", "\n", "mu->setVal(0);\n", "cout << \"is mu=0 in the interval? \" << plInt->IsInInterval(RooArgSet(*mu)) << endl;" ] }, { "cell_type": "markdown", "id": "b147c0d7", "metadata": {}, "source": [ "make a reasonable style" ] }, { "cell_type": "code", "execution_count": null, "id": "fab51952", "metadata": { "collapsed": false }, "outputs": [], "source": [ "gStyle->SetCanvasColor(0);\n", "gStyle->SetCanvasBorderMode(0);\n", "gStyle->SetPadBorderMode(0);\n", "gStyle->SetPadColor(0);\n", "gStyle->SetCanvasColor(0);\n", "gStyle->SetTitleFillColor(0);\n", "gStyle->SetFillColor(0);\n", "gStyle->SetFrameFillColor(0);\n", "gStyle->SetStatColor(0);" ] }, { "cell_type": "markdown", "id": "6ce96a5b", "metadata": {}, "source": [ "some plots" ] }, { "cell_type": "code", "execution_count": null, "id": "8fdc4525", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TCanvas *canvas = new TCanvas(\"canvas\");\n", "canvas->Divide(2, 2);" ] }, { "cell_type": "markdown", "id": "a7b93fc3", "metadata": {}, "source": [ "plot the data" ] }, { "cell_type": "code", "execution_count": null, "id": "d2a477e4", "metadata": { "collapsed": false }, "outputs": [], "source": [ "canvas->cd(1);\n", "RooPlot *frame = x->frame();\n", "data->plotOn(frame);\n", "data->statOn(frame);\n", "frame->Draw();" ] }, { "cell_type": "markdown", "id": "549c7ff0", "metadata": {}, "source": [ "plot the profile likelihood" ] }, { "cell_type": "code", "execution_count": null, "id": "c7e3d4bb", "metadata": { "collapsed": false }, "outputs": [], "source": [ "canvas->cd(2);\n", "LikelihoodIntervalPlot plot(plInt);\n", "plot.Draw();" ] }, { "cell_type": "markdown", "id": "8ddbc98c", "metadata": {}, "source": [ "plot the MCMC interval" ] }, { "cell_type": "code", "execution_count": null, "id": "3ab0d95a", "metadata": { "collapsed": false }, "outputs": [], "source": [ "canvas->cd(3);\n", "MCMCIntervalPlot *mcPlot = new MCMCIntervalPlot(*mcInt);\n", "mcPlot->SetLineColor(kGreen);\n", "mcPlot->SetLineWidth(2);\n", "mcPlot->Draw();\n", "\n", "canvas->cd(4);\n", "RooPlot *bcPlot = bc.GetPosteriorPlot();\n", "bcPlot->Draw();\n", "\n", "canvas->Update();\n", "\n", "t.Stop();\n", "t.Print();" ] }, { "cell_type": "markdown", "id": "dcfc9a35", "metadata": {}, "source": [ "Draw all canvases " ] }, { "cell_type": "code", "execution_count": null, "id": "b3c40010", "metadata": { "collapsed": false }, "outputs": [], "source": [ "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 }