Logo ROOT   6.07/09
Reference Guide
rf610_visualerror.C File Reference
Tutorials » RooFit Tutorials

Detailed Description

View in nbviewer Open in SWAN 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #610

Visualization of errors from a covariance matrix

pict1_rf610_visualerror.C.png
Processing /mnt/vdb/lsf/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roofit/rf610_visualerror.C...
RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (sig,bkg)
**********
** 1 **SET PRINT 1
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 fsig 3.30000e-01 1.00000e-01 0.00000e+00 1.00000e+00
2 m 0.00000e+00 2.00000e+00 -1.00000e+01 1.00000e+01
3 m2 -1.00000e+00 2.00000e+00 -1.00000e+01 1.00000e+01
4 s 2.00000e+00 5.00000e-01 1.00000e+00 5.00000e+01
5 s2 6.00000e+00 2.50000e+00 1.00000e+00 5.00000e+01
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 1
**********
**********
** 5 **SET STR 1
**********
NOW USING STRATEGY 1: TRY TO BALANCE SPEED AGAINST RELIABILITY
**********
** 6 **MIGRAD 2500 1
**********
FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
START MIGRAD MINIMIZATION. STRATEGY 1. CONVERGENCE WHEN EDM .LT. 1.00e-03
FCN=2770.05 FROM MIGRAD STATUS=INITIATE 20 CALLS 21 TOTAL
EDM= unknown STRATEGY= 1 NO ERROR MATRIX
EXT PARAMETER CURRENT GUESS STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 fsig 3.30000e-01 1.00000e-01 2.14988e-01 -9.80006e+00
2 m 0.00000e+00 2.00000e+00 2.01358e-01 -3.95919e+01
3 m2 -1.00000e+00 2.00000e+00 2.02430e-01 3.98515e+01
4 s 2.00000e+00 5.00000e-01 7.46809e-02 2.95415e+01
5 s2 6.00000e+00 2.50000e+00 1.74125e-01 4.56667e+01
ERR DEF= 0.5
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=2767.65 FROM MIGRAD STATUS=CONVERGED 125 CALLS 126 TOTAL
EDM=1.56678e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 fsig 2.98883e-01 6.74303e-02 2.39863e-03 7.00359e-03
2 m 3.08816e-01 2.09098e-01 6.83546e-04 -2.93512e-02
3 m2 -1.31219e+00 3.64277e-01 9.46152e-04 8.97852e-02
4 s 1.78229e+00 2.51997e-01 9.25951e-04 -2.49363e-02
5 s2 5.51197e+00 4.85498e-01 6.94615e-04 1.27501e-01
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 5 ERR DEF=0.5
4.580e-03 -3.800e-03 -1.557e-02 1.331e-02 2.642e-02
-3.800e-03 4.373e-02 -3.141e-03 -1.195e-02 -2.959e-02
-1.557e-02 -3.141e-03 1.328e-01 -4.509e-02 -1.100e-01
1.331e-02 -1.195e-02 -4.509e-02 6.354e-02 7.253e-02
2.642e-02 -2.959e-02 -1.100e-01 7.253e-02 2.358e-01
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4 5
1 0.89430 1.000 -0.269 -0.632 0.780 0.804
2 0.43384 -0.269 1.000 -0.041 -0.227 -0.291
3 0.70478 -0.632 -0.041 1.000 -0.491 -0.621
4 0.78303 0.780 -0.227 -0.491 1.000 0.593
5 0.82883 0.804 -0.291 -0.621 0.593 1.000
**********
** 7 **SET ERR 0.5
**********
**********
** 8 **SET PRINT 1
**********
**********
** 9 **HESSE 2500
**********
COVARIANCE MATRIX CALCULATED SUCCESSFULLY
FCN=2767.65 FROM HESSE STATUS=OK 31 CALLS 157 TOTAL
EDM=1.56685e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER INTERNAL INTERNAL
NO. NAME VALUE ERROR STEP SIZE VALUE
1 fsig 2.98883e-01 6.78952e-02 4.79727e-04 -4.13956e-01
2 m 3.08816e-01 2.09026e-01 1.36709e-04 3.08865e-02
3 m2 -1.31219e+00 3.67042e-01 1.89230e-04 -1.31599e-01
4 s 1.78229e+00 2.53102e-01 1.85190e-04 -1.31741e+00
5 s2 5.51197e+00 4.89300e-01 1.38923e-04 -9.54177e-01
ERR DEF= 0.5
EXTERNAL ERROR MATRIX. NDIM= 25 NPAR= 5 ERR DEF=0.5
4.644e-03 -3.811e-03 -1.596e-02 1.350e-02 2.692e-02
-3.811e-03 4.370e-02 -2.970e-03 -1.206e-02 -2.967e-02
-1.596e-02 -2.970e-03 1.348e-01 -4.619e-02 -1.129e-01
1.350e-02 -1.206e-02 -4.619e-02 6.410e-02 7.402e-02
2.692e-02 -2.967e-02 -1.129e-01 7.402e-02 2.395e-01
PARAMETER CORRELATION COEFFICIENTS
NO. GLOBAL 1 2 3 4 5
1 0.89584 1.000 -0.268 -0.638 0.782 0.807
2 0.43319 -0.268 1.000 -0.039 -0.228 -0.290
3 0.71012 -0.638 -0.039 1.000 -0.497 -0.628
4 0.78518 0.782 -0.228 -0.497 1.000 0.597
5 0.83175 0.807 -0.290 -0.628 0.597 1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) INFO: visualizing 1-sigma uncertainties in parameters (m,s,fsig,m2,s2) from fit result fitresult_model_genData using 315 samplings.
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) INFO: visualizing 1-sigma uncertainties in parameters (m,s,fsig,m2,s2) from fit result fitresult_model_genData using 315 samplings.
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooAddPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TAxis.h"
using namespace RooFit ;
void rf610_visualerror()
{
// S e t u p e x a m p l e f i t
// ---------------------------------------
// Create sum of two Gaussians p.d.f. with factory
RooRealVar x("x","x",-10,10) ;
RooRealVar m("m","m",0,-10,10) ;
RooRealVar s("s","s",2,1,50) ;
RooGaussian sig("sig","sig",x,m,s) ;
RooRealVar m2("m2","m2",-1,-10,10) ;
RooRealVar s2("s2","s2",6,1,50) ;
RooGaussian bkg("bkg","bkg",x,m2,s2) ;
RooRealVar fsig("fsig","fsig",0.33,0,1) ;
RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;
// Create binned dataset
x.setBins(25) ;
RooAbsData* d = model.generateBinned(x,1000) ;
// Perform fit and save fit result
RooFitResult* r = model.fitTo(*d,Save()) ;
// V i s u a l i z e f i t e r r o r
// -------------------------------------
// Make plot frame
RooPlot* frame = x.frame(Bins(40),Title("P.d.f with visualized 1-sigma error band")) ;
d->plotOn(frame) ;
// Visualize 1-sigma error encoded in fit result 'r' as orange band using linear error propagation
// 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, where 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
//
model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange)) ;
// Calculate error using sampling method and visualize as dashed red line.
//
// In this method a number of curves is calculated with variations of the parameter values, as sampled
// from a multi-variate Gaussian p.d.f. 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 value of x, where 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, and 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, but may take (much) longer to calculate
model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed)) ;
// 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, efficiencies etc..)
model.plotOn(frame,VisualizeError(*r,1),FillColor(kOrange),Components("bkg")) ;
model.plotOn(frame,VisualizeError(*r,1,kFALSE),DrawOption("L"),LineWidth(2),LineColor(kRed),Components("bkg"),LineStyle(kDashed)) ;
// Overlay central value
model.plotOn(frame) ;
model.plotOn(frame,Components("bkg"),LineStyle(kDashed)) ;
d->plotOn(frame) ;
frame->SetMinimum(0) ;
// V i s u a l i z e p a r t i a l f i t e r r o r
// ------------------------------------------------------
// Make plot frame
RooPlot* 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,V22 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'), and V22bar
// is the Shur complement of V22, calculated 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
model.plotOn(frame2,VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;
model.plotOn(frame2,Components("bkg"),VisualizeError(*r,RooArgSet(m,m2),2),FillColor(kCyan)) ;
model.plotOn(frame2) ;
model.plotOn(frame2,Components("bkg"),LineStyle(kDashed)) ;
frame2->SetMinimum(0) ;
// Make plot frame
RooPlot* 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
model.plotOn(frame3,VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;
model.plotOn(frame3,Components("bkg"),VisualizeError(*r,RooArgSet(s,s2),2),FillColor(kGreen)) ;
model.plotOn(frame3) ;
model.plotOn(frame3,Components("bkg"),LineStyle(kDashed)) ;
frame3->SetMinimum(0) ;
// Make plot frame
RooPlot* 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
model.plotOn(frame4,VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;
model.plotOn(frame4,Components("bkg"),VisualizeError(*r,RooArgSet(fsig),2),FillColor(kMagenta)) ;
model.plotOn(frame4) ;
model.plotOn(frame4,Components("bkg"),LineStyle(kDashed)) ;
frame4->SetMinimum(0) ;
TCanvas* c = new TCanvas("rf610_visualerror","rf610_visualerror",800,800) ;
c->Divide(2,2) ;
c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;
c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;
c->cd(4) ; gPad->SetLeftMargin(0.15) ; frame4->GetYaxis()->SetTitleOffset(1.6) ; frame4->Draw() ;
}
Author
04/2009 - Wouter Verkerke

Definition in file rf610_visualerror.C.