'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #601
Interactive minimization with MINUIT
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 Tuesday, March 19, 2024 at 07:16 PM.
import ROOT
Observable
x = ROOT.RooRealVar("x", "x", -20, 20)
Model (intentional strong correlations)
mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0)
sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 3)
g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)
sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 6.0)
g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)
frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
model = ROOT.RooAddPdf("model", "model", [g1, g2], [frac])
[#0] WARNING:InputArguments -- The parameter 'sigma_g1' with range [-inf, inf] of the RooGaussian 'g1' exceeds the safe range of (0, inf). Advise to limit its range.
Generate 1000 events
data = model.generate({x}, 1000)
Construct unbinned likelihood of model w.r.t. data
nll = model.createNLL(data)
[#1] INFO:Fitting -- RooAbsPdf::fitTo(model) fixing normalization set for coefficient determination to observables in data [#1] INFO:Fitting -- using CPU computation library compiled with -mavx2
Create MINUIT interface object
m = ROOT.RooMinimizer(nll)
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_modelData) Summation contains a RooNLLVar, using its error level
Activate verbose logging of MINUIT parameter space stepping
m.setVerbose(True)
Call MIGRAD to minimize the likelihood
m.migrad()
0
[#0] WARNING:Minimization -- RooAbsMinimizerFcn::synchronize: WARNING: no initial error estimate available for frac: using 0.1 [#0] WARNING:Minimization -- RooAbsMinimizerFcn::synchronize: WARNING: no initial error estimate available for sigma_g2: using 0.3 Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1 prevFCN = 2660.220684 frac=0.5036, prevFCN = 2660.181264 frac=0.4964, prevFCN = 2660.261875 frac=0.5, sigma_g2=4.011, prevFCN = 2660.278974 sigma_g2=3.989, prevFCN = 2660.167705 sigma_g2=4.005, prevFCN = 2660.248509 sigma_g2=3.995, prevFCN = 2660.194127 frac=0.5812, sigma_g2=3.889, prevFCN = 2660.146969 frac=0.5429, sigma_g2=3.941, prevFCN = 2659.83839 frac=0.5459, prevFCN = 2659.836693 frac=0.5398, prevFCN = 2659.841351 frac=0.5429, sigma_g2=3.946, prevFCN = 2659.835035 sigma_g2=3.936, prevFCN = 2659.842919 frac=0.5497, sigma_g2=3.955, prevFCN = 2659.823248 frac=0.5767, sigma_g2=4.011, prevFCN = 2659.774616 frac=0.6314, sigma_g2=4.128, prevFCN = 2659.73914 frac=0.6266, sigma_g2=4.117, prevFCN = 2659.738319 frac=0.6296, prevFCN = 2659.740343 frac=0.6237, prevFCN = 2659.737969 frac=0.6266, sigma_g2=4.123, prevFCN = 2659.737996 sigma_g2=4.112, prevFCN = 2659.739643 frac=0.6227, sigma_g2=4.114, prevFCN = 2659.737959 frac=0.6236, sigma_g2=4.115, prevFCN = 2659.737923 frac=0.6262, prevFCN = 2659.738617 frac=0.621, prevFCN = 2659.738491 frac=0.6236, sigma_g2=4.121, prevFCN = 2659.738401 sigma_g2=4.108, prevFCN = 2659.738723 sigma_g2=4.115, prevFCN = 2659.737923 frac=0.6262, prevFCN = 2659.738617 frac=0.621, prevFCN = 2659.738491 frac=0.6236, sigma_g2=4.121, prevFCN = 2659.738401 sigma_g2=4.108, prevFCN = 2659.738723 frac=0.6241, sigma_g2=4.115, prevFCN = 2659.737961 frac=0.6231, prevFCN = 2659.737935 frac=0.6236, sigma_g2=4.116, prevFCN = 2659.737916 sigma_g2=4.113, prevFCN = 2659.737981 frac=0.6262, sigma_g2=4.121, prevFCN = 2659.737949 Minuit2Minimizer : Valid minimum - status = 0 FVAL = 2659.73792283928833 Edm = 2.4027910430005747e-05 Nfcn = 37 frac = 0.6236 +/- 0.164 (limited) sigma_g2 = 4.115 +/- 0.405 (limited) frac=0.6236, sigma_g2=4.115,
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.220684 Edm = 0.7499658139 NCalls = 7 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.220684 Edm : 0.7499658139 Internal parameters: [ 0 -0.3398369095] Internal gradient : [ -5.61967122 7.285345928] Internal covariance matrix: [[ 0.058086658 0] [ 0 0.021957944]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.220684 Edm = 0.7499658139 NCalls = 7 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2659.83839 Edm = 0.008262338757 NCalls = 13 Info in <Minuit2>: VariableMetricBuilder 2 - FCN = 2659.738319 Edm = 0.0004391718851 NCalls = 21 Info in <Minuit2>: VariableMetricBuilder 3 - FCN = 2659.737923 Edm = 2.063673982e-05 NCalls = 27 Info in <Minuit2>: VariableMetricBuilder After Hessian Info in <Minuit2>: VariableMetricBuilder 4 - FCN = 2659.737923 Edm = 2.402791043e-05 NCalls = 37
Print values of all parameters, reflect values (and error estimates) that are back propagated from MINUIT
model.getParameters({x}).Print("s")
1) RooRealVar:: frac = 0.6236 +/- 0.164 2) RooRealVar:: mean = 0 3) RooRealVar:: sigma_g1 = 3 4) RooRealVar:: sigma_g2 = 4.115 +/- 0.405
Disable verbose logging
m.setVerbose(False)
Run HESSE to calculate errors from d2L/dp2
m.hesse()
0
Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 1000 Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate
Print value (and error) of sigma_g2 parameter, reflects value and error back propagated from MINUIT
sigma_g2.Print()
RooRealVar::sigma_g2 = 4.115 +/- 0.4057 L(3 - 6)
Run MINOS on sigma_g2 parameter only
m.minos({sigma_g2})
0
****************************************************************************************************** Minuit2Minimizer::GetMinosError - Run MINOS LOWER error for parameter #1 : sigma_g2 using max-calls 1000, tolerance 1 ****************************************************************************************************** Minuit2Minimizer::GetMinosError - Run MINOS UPPER error for parameter #1 : sigma_g2 using max-calls 1000, tolerance 1 Minos: Lower error for parameter sigma_g2 : -0.3794 Minos: Upper error for parameter sigma_g2 : 0.4574
Info in <Minuit2>: MnMinos Determination of lower Minos error for parameter 1 Info in <Minuit2>: MnFunctionCross Run Migrad with fixed parameters: Pos 1: sigma_g2 = 3.70891 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.55118 Edm = 0.113308465 NCalls = 3 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.55118 Edm : 0.113308465 Internal parameters: [ -0.06226626228] Internal gradient : [ 3.194420786] Internal covariance matrix: [[ 0.044415863]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.55118 Edm = 0.113308465 NCalls = 3 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.315058 Edm = 0.0002150022484 NCalls = 7 Info in <Minuit2>: MnFunctionCross Result after Migrad FCN = 2660.315058 Edm = 0.0002150022484 NCalls = 7 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.3972435494 | 0.1056591218 1 | sigma_g2 | fixed | 3.70891435 | 0.4056662034 Info in <Minuit2>: MnFunctionCross Run Migrad again (2nd) with fixed parameters: Pos 1: sigma_g2 = 3.737 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.250576 Edm = 0.01729034995 NCalls = 3 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.250576 Edm : 0.01729034995 Internal parameters: [ -0.2069877666] Internal gradient : [ -0.8544275172] Internal covariance matrix: [[ 0.094735624]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.250576 Edm = 0.01729034995 NCalls = 3 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.232876 Edm = 4.007072232e-06 NCalls = 6 Info in <Minuit2>: MnFunctionCross Result after 2nd Migrad FCN = 2660.232876 Edm = 4.007072232e-06 NCalls = 6 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.4171265177 | 0.1072770084 1 | sigma_g2 | fixed | 3.736995131 | 0.4056662034 Info in <Minuit2>: MnFunctionCross Run Migrad again (3rd) with fixed parameters: Pos 1: sigma_g2 = 3.73527 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.237741 Edm = 5.406683871e-05 NCalls = 3 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.237741 Edm : 5.406683871e-05 Internal parameters: [ -0.1665154044] Internal gradient : [ 0.04741411879] Internal covariance matrix: [[ 0.096200117]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.237741 Edm = 5.406683871e-05 NCalls = 3 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.237689 Edm = 2.189992327e-08 NCalls = 6 Info in <Minuit2>: MnFunctionCross Result after Migrad (3rd): FCN = 2660.237689 Edm = 2.189992327e-08 NCalls = 6 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.4160021957 | 0.1061893687 1 | sigma_g2 | fixed | 3.735270747 | 0.4056662034 Info in <Minuit2>: MnMinos end of Minos scan for low interval for parameter sigma_g2 Info in <Minuit2>: MnMinos Determination of upper Minos error for parameter 1 Info in <Minuit2>: MnFunctionCross Run Migrad with fixed parameters: Pos 1: sigma_g2 = 4.52025 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.321661 Edm = 0.2677306116 NCalls = 5 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.321661 Edm : 0.2677306116 Internal parameters: [ 0.5618487476] Internal gradient : [ 4.910321861] Internal covariance matrix: [[ 0.044415863]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.321661 Edm = 0.2677306116 NCalls = 5 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.139576 Edm = 5.875658783e-08 NCalls = 9 Info in <Minuit2>: MnFunctionCross Result after Migrad FCN = 2660.139576 Edm = 5.875658783e-08 NCalls = 9 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.7343787943 | 0.0540577002 1 | sigma_g2 | fixed | 4.520246757 | 0.4056662034 Info in <Minuit2>: MnFunctionCross Run Migrad again (2nd) with fixed parameters: Pos 1: sigma_g2 = 4.56719 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.244047 Edm = 0.01621250066 NCalls = 3 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.244047 Edm : 0.01621250066 Internal parameters: [ 0.4878837384] Internal gradient : [ -1.467518302] Internal covariance matrix: [[ 0.030112232]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.244047 Edm = 0.01621250066 NCalls = 3 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.228567 Edm = 3.216983647e-05 NCalls = 6 Info in <Minuit2>: MnFunctionCross Result after 2nd Migrad FCN = 2660.228567 Edm = 3.216983647e-05 NCalls = 6 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.7440794003 | 0.05226504471 1 | sigma_g2 | fixed | 4.567194795 | 0.4056662034 Info in <Minuit2>: MnFunctionCross Run Migrad again (3rd) with fixed parameters: Pos 1: sigma_g2 = 4.57213 Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2660.238339 Edm = 5.179530532e-05 NCalls = 3 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2660.238339 Edm : 5.179530532e-05 Internal parameters: [ 0.5099788641] Internal gradient : [ -0.08477490437] Internal covariance matrix: [[ 0.028828084]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.0005 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2660.238339 Edm = 5.179530532e-05 NCalls = 3 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2660.238288 Edm = 1.944203473e-10 NCalls = 6 Info in <Minuit2>: MnFunctionCross Result after Migrad (3rd): FCN = 2660.238288 Edm = 1.944203473e-10 NCalls = 6 Pos | Name | type | Value | Error +/- 0 | frac | limited | 0.7446124489 | 0.05217900467 1 | sigma_g2 | fixed | 4.572130423 | 0.4056662034 Info in <Minuit2>: MnMinos end of Minos scan for up interval for parameter sigma_g2
Print value (and error) of sigma_g2 parameter, reflects value and error back propagated from MINUIT
sigma_g2.Print()
RooRealVar::sigma_g2 = 4.115 +/- (-0.3794,0.4574) L(3 - 6)
Save a snapshot of the fit result. ROOT.This object contains the initial fit parameters, final fit parameters, complete correlation matrix, EDM, minimized FCN , last MINUIT status code and the number of times the ROOT.RooFit function object has indicated evaluation problems (e.g. zero probabilities during likelihood evaluation)
r = m.save()
Make contour plot of mx vs sx at 1,2, sigma
frame = m.contour(frac, sigma_g2, 1, 2, 3)
frame.SetTitle("Contour plot")
Info in <Minuit2>: Minuit2Minimizer::Contour Computing contours - 0.5 Info in <Minuit2>: Minuit2Minimizer::Contour Computing contours - 2 Info in <Minuit2>: Minuit2Minimizer::Contour Computing contours - 4.5
Print the fit result snapshot
r.Print("v")
RooFitResult: minimized FCN value: 2660, estimated distance to minimum: 2.409e-05 covariance matrix quality: Full, accurate covariance matrix Status : MIGRAD=0 HESSE=0 MINOS=0 Constant Parameter Value -------------------- ------------ mean 0.0000e+00 sigma_g1 3.0000e+00 Floating Parameter InitialValue FinalValue (+HiError,-LoError) GblCorr. -------------------- ------------ ---------------------------------- -------- frac 5.0000e-01 6.2360e-01 +/- 1.64e-01 <none> sigma_g2 4.0000e+00 4.1146e+00 (+4.57e-01,-3.79e-01) <none>
At any moment you can manually change the value of a (constant) parameter
mean.setVal(0.3)
Rerun MIGRAD,HESSE
m.migrad()
m.hesse()
frac.Print()
Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1 Minuit2Minimizer : Valid minimum - status = 0 FVAL = 2663.35774508719987 Edm = 9.56369873007578491e-05 Nfcn = 38 frac = 0.5655 +/- 0.1961 (limited) sigma_g2 = 4.005 +/- 0.3917 (limited) RooRealVar::frac = 0.5655 +/- 0.1961 L(0 - 1)
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2663.411598 Edm = 0.0145493914 NCalls = 9 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2663.411598 Edm : 0.0145493914 Internal parameters: [ 0.2497912427 -0.2598610789] Internal gradient : [ 1.120322622 -0.2622810121] Internal covariance matrix: [[ 0.044738711 0] [ 0 0.029727437]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2663.411598 Edm = 0.0145493914 NCalls = 9 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2663.400666 Edm = 0.005796742296 NCalls = 15 Info in <Minuit2>: VariableMetricBuilder 2 - FCN = 2663.35958 Edm = 0.00172901389 NCalls = 22 Info in <Minuit2>: VariableMetricBuilder 3 - FCN = 2663.357745 Edm = 8.583255993e-05 NCalls = 28 Info in <Minuit2>: VariableMetricBuilder After Hessian Info in <Minuit2>: VariableMetricBuilder 4 - FCN = 2663.357745 Edm = 9.56369873e-05 NCalls = 38 Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 1000 Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate
Now fix sigma_g2
sigma_g2.setConstant(True)
Rerun MIGRAD,HESSE
m.migrad()
m.hesse()
frac.Print()
c = ROOT.TCanvas("rf601_intminuit", "rf601_intminuit", 600, 600)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.4)
frame.Draw()
c.SaveAs("rf601_intminuit.png")
Minuit2Minimizer: Minimize with max-calls 1000 convergence for edm < 1 strategy 1 Minuit2Minimizer : Valid minimum - status = 0 FVAL = 2663.35774092453448 Edm = 1.35713798793672428e-10 Nfcn = 15 frac = 0.5652 +/- 0.08029 (limited) sigma_g2 = 4.005 (fixed) RooRealVar::frac = 0.5652 +/- 0.08029 L(0 - 1)
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 2663.357745 Edm = 4.211167051e-06 NCalls = 5 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 2663.357745 Edm : 4.211167051e-06 Internal parameters: [ 0.1313437863] Internal gradient : [ 0.01784069115] Internal covariance matrix: [[ 0.052922349]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.001 with call limit = 1000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 2663.357745 Edm = 4.211167051e-06 NCalls = 5 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 2663.357741 Edm = 6.315677937e-14 NCalls = 8 Info in <Minuit2>: VariableMetricBuilder After Hessian Info in <Minuit2>: VariableMetricBuilder 2 - FCN = 2663.357741 Edm = 1.357137988e-10 NCalls = 15 Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 1000 Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate Info in <TCanvas::Print>: png file rf601_intminuit.png has been created
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()