Logo ROOT  
Reference Guide
rf610_visualerror.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Likelihood and minimization: visualization of errors from a covariance matrix
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import ROOT
12
13
14# Setup example fit
15# ---------------------------------------
16
17# Create sum of two Gaussians pdf with factory
18x = ROOT.RooRealVar("x", "x", -10, 10)
19
20m = ROOT.RooRealVar("m", "m", 0, -10, 10)
21s = ROOT.RooRealVar("s", "s", 2, 1, 50)
22sig = ROOT.RooGaussian("sig", "sig", x, m, s)
23
24m2 = ROOT.RooRealVar("m2", "m2", -1, -10, 10)
25s2 = ROOT.RooRealVar("s2", "s2", 6, 1, 50)
26bkg = ROOT.RooGaussian("bkg", "bkg", x, m2, s2)
27
28fsig = ROOT.RooRealVar("fsig", "fsig", 0.33, 0, 1)
29model = ROOT.RooAddPdf("model", "model", [sig, bkg], [fsig])
30
31# Create binned dataset
32x.setBins(25)
33d = model.generateBinned({x}, 1000)
34
35# Perform fit and save fit result
36r = model.fitTo(d, Save=True)
37
38# Visualize fit error
39# -------------------------------------
40
41# Make plot frame
42frame = x.frame(Bins=40, Title="P.d.f with visualized 1-sigma error band")
43d.plotOn(frame)
44
45# Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
46# ROOT.This results in an error band that is by construction symmetric
47#
48# The linear error is calculated as
49# error(x) = Z* F_a(x) * Corr(a,a') F_a'(x)
50#
51# where F_a(x) = [ f(x,a+da) - f(x,a-da) ] / 2,
52#
53# with f(x) = the plotted curve
54# 'da' = error taken from the fit result
55# Corr(a,a') = the correlation matrix from the fit result
56# Z = requested significance 'Z sigma band'
57#
58# The linear method is fast (required 2*N evaluations of the curve, N is the number of parameters),
59# but may not be accurate in the presence of strong correlations (~>0.9) and at Z>2 due to linear and
60# Gaussian approximations made
61#
62model.plotOn(frame, VisualizeError=(r, 1), FillColor="kOrange")
63
64# Calculate error using sampling method and visualize as dashed red line.
65#
66# In self method a number of curves is calculated with variations of the parameter values, sampled
67# from a multi-variate Gaussian pdf that is constructed from the fit results covariance matrix.
68# The error(x) is determined by calculating a central interval that capture N% of the variations
69# for each valye of x, N% is controlled by Z (i.e. Z=1 gives N=68%). The number of sampling curves
70# is chosen to be such that at least 100 curves are expected to be outside the N% interval, is minimally
71# 100 (e.g. Z=1.Ncurve=356, Z=2.Ncurve=2156)) Intervals from the sampling method can be asymmetric,
72# and may perform better in the presence of strong correlations, may take
73# (much) longer to calculate
74model.plotOn(frame, VisualizeError=(r, 1, False), DrawOption="L", LineWidth=2, LineColor="r")
75
76# Perform the same type of error visualization on the background component only.
77# The VisualizeError() option can generally applied to _any_ kind of
78# plot (components, asymmetries, etc..)
79model.plotOn(frame, VisualizeError=(r, 1), FillColor="kOrange", Components="bkg")
80model.plotOn(
81 frame,
82 VisualizeError=(r, 1, False),
83 DrawOption="L",
84 LineWidth=2,
85 LineColor="r",
86 Components="bkg",
87 LineStyle="--",
88)
89
90# Overlay central value
91model.plotOn(frame)
92model.plotOn(frame, Components="bkg", LineStyle="--")
93d.plotOn(frame)
94frame.SetMinimum(0)
95
96# Visualize partial fit error
97# ------------------------------------------------------
98
99# Make plot frame
100frame2 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from (m,m2)")
101
102# Visualize partial error. For partial error visualization the covariance matrix is first reduced as follows
103# ___ -1
104# Vred = V22 = V11 - V12 * V22 * V21
105#
106# Where V11,V12,V21, represent a block decomposition of the covariance matrix into observables that
107# are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), V22bar
108# is the Shur complement of V22, as shown above
109#
110# (Note that Vred is _not_ a simple sub-matrix of V)
111
112# Propagate partial error due to shape parameters (m,m2) using linear and
113# sampling method
114model.plotOn(frame2, VisualizeError=(r, {m, m2}, 2), FillColor="c")
115model.plotOn(frame2, Components="bkg", VisualizeError=(r, {m, m2}, 2), FillColor="c")
116
117model.plotOn(frame2)
118model.plotOn(frame2, Components="bkg", LineStyle="--")
119frame2.SetMinimum(0)
120
121# Make plot frame
122frame3 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from (s,s2)")
123
124# Propagate partial error due to yield parameter using linear and sampling
125# method
126model.plotOn(frame3, VisualizeError=(r, {s, s2}, 2), FillColor="g")
127model.plotOn(frame3, Components="bkg", VisualizeError=(r, {fsig}, 2), FillColor="g")
128
129model.plotOn(frame3)
130model.plotOn(frame3, Components="bkg", LineStyle="--")
131frame3.SetMinimum(0)
132
133# Make plot frame
134frame4 = x.frame(Bins=40, Title="Visualization of 2-sigma partial error from fsig")
135
136# Propagate partial error due to yield parameter using linear and sampling
137# method
138model.plotOn(frame4, VisualizeError=(r, {fsig}, 2), FillColor="m")
139model.plotOn(frame4, Components="bkg", VisualizeError=(r, {fsig}, 2), FillColor="m")
140
141model.plotOn(frame4)
142model.plotOn(frame4, Components="bkg", LineStyle="--")
143frame4.SetMinimum(0)
144
145c = ROOT.TCanvas("rf610_visualerror", "rf610_visualerror", 800, 800)
146c.Divide(2, 2)
147c.cd(1)
148ROOT.gPad.SetLeftMargin(0.15)
149frame.GetYaxis().SetTitleOffset(1.4)
150frame.Draw()
151c.cd(2)
152ROOT.gPad.SetLeftMargin(0.15)
153frame2.GetYaxis().SetTitleOffset(1.6)
154frame2.Draw()
155c.cd(3)
156ROOT.gPad.SetLeftMargin(0.15)
157frame3.GetYaxis().SetTitleOffset(1.6)
158frame3.Draw()
159c.cd(4)
160ROOT.gPad.SetLeftMargin(0.15)
161frame4.GetYaxis().SetTitleOffset(1.6)
162frame4.Draw()
163
164c.SaveAs("rf610_visualerror.png")