Likelihood and minimization: visualization of errors from a covariance matrix

**Author:** Clemens Lange, Wouter Verkerke (C++ version)

*This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Friday, May 13, 2022 at 09:41 AM.*

In [1]:

```
import ROOT
```

Create sum of two Gaussians pdf with factory

In [2]:

```
x = ROOT.RooRealVar("x", "x", -10, 10)
m = ROOT.RooRealVar("m", "m", 0, -10, 10)
s = ROOT.RooRealVar("s", "s", 2, 1, 50)
sig = ROOT.RooGaussian("sig", "sig", x, m, s)
m2 = ROOT.RooRealVar("m2", "m2", -1, -10, 10)
s2 = ROOT.RooRealVar("s2", "s2", 6, 1, 50)
bkg = ROOT.RooGaussian("bkg", "bkg", x, m2, s2)
fsig = ROOT.RooRealVar("fsig", "fsig", 0.33, 0, 1)
model = ROOT.RooAddPdf("model", "model", [sig, bkg], [fsig])
```

Create binned dataset

In [3]:

```
x.setBins(25)
d = model.generateBinned({x}, 1000)
```

Perform fit and save fit result

In [4]:

```
r = model.fitTo(d, Save=True)
```

Make plot frame

In [5]:

```
frame = x.frame(Bins=40, Title="P.d.f with visualized 1-sigma error band")
d.plotOn(frame)
```

Out[5]:

Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation ROOT.This results in an error band that is by construction symmetric

The linear error is calculated as
error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)

where F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2,

```
with f(x) = the plotted curve
'da' = error taken from the fit result
Corr(a,a') = the correlation matrix from the fit result
Z = requested significance 'Z sigma band'
```

The linear method is fast (required 2*N evaluations of the curve, N is the number of parameters), but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and Gaussian approximations made

In [6]:

```
model.plotOn(frame, VisualizeError=(r, 1), FillColor="kOrange")
```

Out[6]:

Calculate error using sampling method and visualize as dashed red line.

In self method a number of curves is calculated with variations of the parameter values, sampled from a multi-variate Gaussian pdf that is constructed from the fit results covariance matrix. The error(x) is determined by calculating a central interval that capture N% of the variations for each valye of x, N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves is chosen to be such that at least 100 curves are expected to be outside the N% interval, is minimally 100 (e.g. Z=1.Ncurve=356, Z=2.Ncurve=2156)) Intervals from the sampling method can be asymmetric, and may perform better in the presence of strong correlations, may take (much) longer to calculate

In [7]:

```
model.plotOn(frame, VisualizeError=(r, 1, False), DrawOption="L", LineWidth=2, LineColor="r")
```

Out[7]:

Perform the same type of error visualization on the background component only.
The VisualizeError() option can generally applied to *any* kind of
plot (components, asymmetries, etc..)

In [8]:

```
model.plotOn(frame, VisualizeError=(r, 1), FillColor="kOrange", Components="bkg")
model.plotOn(
frame,
VisualizeError=(r, 1, False),
DrawOption="L",
LineWidth=2,
LineColor="r",
Components="bkg",
LineStyle="--",
)
```

Out[8]:

Overlay central value

In [9]:

```
model.plotOn(frame)
model.plotOn(frame, Components="bkg", LineStyle="--")
d.plotOn(frame)
frame.SetMinimum(0)
```

Make plot frame

In [10]:

```
frame2 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from (m,m2)")
```

Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
___ -1
Vred = V22 = V11 - V12 * V22 * V21

Where V11,V12,V21, represent a block decomposition of the covariance matrix into observables that are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), V22bar is the Shur complement of V22, as shown above

(Note that Vred is *not* a simple sub-matrix of V)

Propagate partial error due to shape parameters (m,m2) using linear and sampling method

In [11]:

```
model.plotOn(frame2, VisualizeError=(r, {m, m2}, 2), FillColor="c")
model.plotOn(frame2, Components="bkg", VisualizeError=(r, {m, m2}, 2), FillColor="c")
model.plotOn(frame2)
model.plotOn(frame2, Components="bkg", LineStyle="--")
frame2.SetMinimum(0)
```

Make plot frame

In [12]:

```
frame3 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from (s,s2)")
```

Propagate partial error due to yield parameter using linear and sampling method

In [13]:

```
model.plotOn(frame3, VisualizeError=(r, {s, s2}, 2), FillColor="g")
model.plotOn(frame3, Components="bkg", VisualizeError=(r, {fsig}, 2), FillColor="g")
model.plotOn(frame3)
model.plotOn(frame3, Components="bkg", LineStyle="--")
frame3.SetMinimum(0)
```

Make plot frame

In [14]:

```
frame4 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from fsig")
```

Propagate partial error due to yield parameter using linear and sampling method

In [15]:

```
model.plotOn(frame4, VisualizeError=(r, {fsig}, 2), FillColor="m")
model.plotOn(frame4, Components="bkg", VisualizeError=(r, {fsig}, 2), FillColor="m")
model.plotOn(frame4)
model.plotOn(frame4, Components="bkg", LineStyle="--")
frame4.SetMinimum(0)
c = ROOT.TCanvas("rf610_visualerror", "rf610_visualerror", 800, 800)
c.Divide(2, 2)
c.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.4)
frame.Draw()
c.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
frame2.GetYaxis().SetTitleOffset(1.6)
frame2.Draw()
c.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
frame3.GetYaxis().SetTitleOffset(1.6)
frame3.Draw()
c.cd(4)
ROOT.gPad.SetLeftMargin(0.15)
frame4.GetYaxis().SetTitleOffset(1.6)
frame4.Draw()
c.SaveAs("rf610_visualerror.png")
```

Draw all canvases

In [16]:

```
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()
```