Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
rf601_intminuit.py File Reference

Detailed Description

View in nbviewer Open in SWAN

'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #601

Interactive minimization with MINUIT

import ROOT
# Setup pdf and likelihood
# -----------------------------------------------
# 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])
# Generate 1000 events
data = model.generate({x}, 1000)
# Construct unbinned likelihood of model w.r.t. data
nll = model.createNLL(data)
# Interactive minimization, error analysis
# -------------------------------------------------------------------------------
# Create MINUIT interface object
m = ROOT.RooMinimizer(nll)
# Activate verbose logging of MINUIT parameter space stepping
m.setVerbose(True)
# Call MIGRAD to minimize the likelihood
m.migrad()
# Print values of all parameters, reflect values (and error estimates)
# that are back propagated from MINUIT
model.getParameters({x}).Print("s")
# Disable verbose logging
m.setVerbose(False)
# Run HESSE to calculate errors from d2L/dp2
m.hesse()
# Print value (and error) of sigma_g2 parameter, reflects
# value and error back propagated from MINUIT
sigma_g2.Print()
# Run MINOS on sigma_g2 parameter only
m.minos({sigma_g2})
# Print value (and error) of sigma_g2 parameter, reflects
# value and error back propagated from MINUIT
sigma_g2.Print()
# Saving results, contour plots
# ---------------------------------------------------------
# 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")
# Print the fit result snapshot
r.Print("v")
# Change parameter values, plotting
# -----------------------------------------------------------------
# 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()
# 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")
void Print(GNN_Data &d, std::string txt="")
[#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.
[#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
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_model_modelData) Summation contains a RooNLLVar, using its error level
[#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, 1) RooRealVar:: frac = 0.6236 +/- 0.164
2) RooRealVar:: mean = 0
3) RooRealVar:: sigma_g1 = 3
4) RooRealVar:: sigma_g2 = 4.115 +/- 0.405
RooRealVar::sigma_g2 = 4.115 +/- 0.4057 L(3 - 6)
******************************************************************************************************
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
RooRealVar::sigma_g2 = 4.115 +/- (-0.3794,0.4574) L(3 - 6)
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>
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)
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)
Date
February 2018
Authors
Clemens Lange, Wouter Verkerke (C version)

Definition in file rf601_intminuit.py.