Introduction to Minuit2.jl
This tutorial is based in the iminuit introduction tutorial It is a simple example that shows how to interact with the Minuit2 classes.
You can also download this example as a Jupyter notebook and a plain Julia source file.
Table of contents
- Loading the necessary Julia modules
- Define some data to fit
- Define a cost function
- Create a Minuit object
- Minimize the cost function
- Fit of model with flexible number of parameters
- Parameter uncertainties
- Call asymmetric uncertainty intervals ("Minos errors")
- Built-in plotting
- Summary
Loading the necessary Julia modules
Load the Minuit2 module. We will also use the Polynomials and Plots modules to define cost functions display results.
using Minuit2
using Polynomials
using FHist
using PlotsDefine some data to fit
We generate some data points that follow a linear model with some noise.
# let's make a line model
line(x, a, b) = a + x * b
a_true = 1.0
b_true = 2.0
# let's make some data
x = range(0, 1., 10)
# precomputed random numbers from standard normal distribution
z = [-0.49783783, -0.33041722, -1.71800806, 1.60229399,
1.36682387, -1.15424221, -0.91425267, -0.03395604,
-1.27611719, -0.7004073 ]
sigma_y = 0.1
y = line.(x, a_true, b_true) + sigma_y * z
# Plot with error bars
plot(x, y, yerr=sigma_y, seriestype=:scatter)- Want to estimate parameters (𝑎,𝑏) of line model from data
- Need a
score` which is minimal when model best agrees with data- Sum of residuals squared (least-squares method)
- Negated sum of log-likelihood values (maximum-likelihood method)
- MINUIT always minimizes; negate score function to maximize
- Use
Minuit2to numerically minimize score as function of model parameters
Define a cost function
We will use the least squares method to define the cost function. The cost function is the sum of the residuals squared.
LSQ(a, b) = sum((y - line.(x, a, b)) .^ 2 ./ sigma_y .^ 2)LSQ (generic function with 1 method)Create a Minuit object
Create a Minuit object with the cost function and initial parameter values. We fix the parameter a to an arbitrary value on this first attempt.
m = Minuit(LSQ; a=2, b=5, fix_a=true,
error_a=0.1, error_b=0.1,
errordef=1)Minuit(FCN = LSQ(a, b), X0 = [2, 5], Method = migrad)Minimize the cost function
Minimize the cost function using the migrad method
migrad!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 305.114 │ migrad │ 13 │ 3 │ 229.214 ms │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 1.31072e-14 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬──────┬──────────┬─────────────┬────────┬────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼──────┼──────────┼─────────────┼────────┼────────┼────────┼────────┼───────┼───────┤
│ 1 │ a │ 2.0 │ 0.1 │ │ │ │ │ true │ │
│ 2 │ b │ 0.511055 │ 0.0533114 │ │ │ │ │ │ │
└───┴──────┴──────────┴─────────────┴────────┴────────┴────────┴────────┴───────┴───────┘
┌───┬─────┬─────┐
│ │ a │ b │
├───┼─────┼─────┤
│ a │ 1.0 │ 0.0 │
│ b │ 0.0 │ 1.0 │
└───┴─────┴─────┘
Plot the data and the fitted line
# get parameter values
a_fit, b_fit = m.values
# Plot with error bars
plot(x, y, yerr=sigma_y, seriestype=:scatter)
plot!(x, line.(x, a_fit, b_fit))Lets change the initial value of a and b and run migrad again
m = Minuit(LSQ; a=5, b=5, error_a=0.1, error_b=0.1,
limit_a=(0, Inf), limit_b=(0, 10.), errordef=1)
migrad!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 10.387 │ migrad │ 64 │ 8 │ 426.405 µs │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 1.93139e-6 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬──────┬──────────┬─────────────┬────────┬────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼──────┼──────────┼─────────────┼────────┼────────┼────────┼────────┼───────┼───────┤
│ 1 │ a │ 0.990894 │ 0.0587697 │ │ │ 0.0 │ │ │ │
│ 2 │ b │ 1.94501 │ 0.0990814 │ │ │ 0.0 │ 10.0 │ │ │
└───┴──────┴──────────┴─────────────┴────────┴────────┴────────┴────────┴───────┴───────┘
┌───┬───────────┬───────────┐
│ │ a │ b │
├───┼───────────┼───────────┤
│ a │ 1.0 │ -0.842944 │
│ b │ -0.842944 │ 1.0 │
└───┴───────────┴───────────┘
Plot the data and the fitted line
# get parameter values
a_fit, b_fit = m.values
# plot with error bars
plot(x, y, yerr=sigma_y, seriestype=:scatter)
plot!(x, line.(x, a_fit, b_fit))Fit of model with flexible number of parameters
- Sometimes the model has large or variable number of parameters
- Example: fit a polynomial of degree 2, 3, 4, ... ?
- Minuit2 has alternative interface which passes parameters as
AbstractVectorto the score function
Define a polynomial model
function LSQ_v(par) # par is a vector of parameters
pol = Polynomial(par) # for len(par) == 2 this is a line
sum((y - pol.(x)) .^ 2 ./ sigma_y .^ 2)
end
# This is the order of coefficients in the polynomial (reverse order in np.polyval)
Polynomial([1,2,3,4])Create a Minuit object with the cost function and initial parameter values
m = Minuit(LSQ_v, [5, 5], error=[0.1, 0.1], errordef=1)
migrad!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 10.387 │ migrad │ 32 │ 4 │ 56.65 ms │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 2.20208e-22 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬────────┬──────────┬─────────────┬────────┬────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼────────┼──────────┼─────────────┼────────┼────────┼────────┼────────┼───────┼───────┤
│ 1 │ par[1] │ 0.990966 │ 0.0587754 │ │ │ │ │ │ │
│ 2 │ par[2] │ 1.94494 │ 0.0990867 │ │ │ │ │ │ │
└───┴────────┴──────────┴─────────────┴────────┴────────┴────────┴────────┴───────┴───────┘
┌────────┬───────────┬───────────┐
│ │ par[1] │ par[2] │
├────────┼───────────┼───────────┤
│ par[1] │ 1.0 │ -0.842927 │
│ par[2] │ -0.842927 │ 1.0 │
└────────┴───────────┴───────────┘
names are automatically generated for the parameters, or can be explicitly set
m = Minuit(LSQ_v, [2, 1, 3, 5], error=0.1,
names=("a", "b", "c", "d"), errordef=1)
migrad!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 9.033 │ migrad │ 96 │ 7 │ 481.346 µs │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 1.39033e-19 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬──────┬──────────┬─────────────┬────────┬────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼──────┼──────────┼─────────────┼────────┼────────┼────────┼────────┼───────┼───────┤
│ 1 │ a │ 0.911937 │ 0.0907623 │ │ │ │ │ │ │
│ 2 │ b │ 2.73544 │ 0.831502 │ │ │ │ │ │ │
│ 3 │ c │ -1.50305 │ 1.99888 │ │ │ │ │ │ │
│ 4 │ d │ 0.765475 │ 1.3117 │ │ │ │ │ │ │
└───┴──────┴──────────┴─────────────┴────────┴────────┴────────┴────────┴───────┴───────┘
┌───┬───────────┬───────────┬───────────┬───────────┐
│ │ a │ b │ c │ d │
├───┼───────────┼───────────┼───────────┼───────────┤
│ a │ 1.0 │ -0.757506 │ 0.593218 │ -0.499578 │
│ b │ -0.757506 │ 1.0 │ -0.958508 │ 0.897816 │
│ c │ 0.593218 │ -0.958508 │ 1.0 │ -0.984327 │
│ d │ -0.499578 │ 0.897816 │ -0.984327 │ 1.0 │
└───┴───────────┴───────────┴───────────┴───────────┘
Lets plot the data and the fitted polynomial
# get parameter values
par_fit = m.values
pol = Polynomial(par_fit)
# Plot with error bars
plot(x, y, yerr=sigma_y, seriestype=:scatter)
plot!(x, line.(x, a_fit, b_fit), label="pol2")
plot!(x, pol.(x), label="pol4")Lets check reduced chi2, goodness-of-fit estimate, should be around 1
m.fval / (length(y) - length(par_fit))1.5055320143296262Parameter uncertainties
- Minuit2 can compute symmetric uncertainty intervals ("Hesse errors")
- automatically done during standard minimisation
- to make sure you get accurate errors, call hesse!(m) explicitly after migrad!(m)
- slow, computation time scales with 𝑁par^2
- Minuit2 can also compute asymmetric uncertainty intervals ("Minos errors")
- need to explicitly call m.minos()
- very slow, computation time scales with 𝑁par^2
Call hesse to get parameter errors
hesse!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 9.033 │ migrad │ 119 │ 8 │ 481.346 µs │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 7.48753e-19 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬──────┬──────────┬─────────────┬────────┬────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼──────┼──────────┼─────────────┼────────┼────────┼────────┼────────┼───────┼───────┤
│ 1 │ a │ 0.911937 │ 0.0907645 │ │ │ │ │ │ │
│ 2 │ b │ 2.73544 │ 0.831558 │ │ │ │ │ │ │
│ 3 │ c │ -1.50305 │ 1.99903 │ │ │ │ │ │ │
│ 4 │ d │ 0.765475 │ 1.31179 │ │ │ │ │ │ │
└───┴──────┴──────────┴─────────────┴────────┴────────┴────────┴────────┴───────┴───────┘
┌───┬───────────┬───────────┬───────────┬───────────┐
│ │ a │ b │ c │ d │
├───┼───────────┼───────────┼───────────┼───────────┤
│ a │ 1.0 │ -0.757518 │ 0.593245 │ -0.499613 │
│ b │ -0.757518 │ 1.0 │ -0.958513 │ 0.89783 │
│ c │ 0.593245 │ -0.958513 │ 1.0 │ -0.984329 │
│ d │ -0.499613 │ 0.89783 │ -0.984329 │ 1.0 │
└───┴───────────┴───────────┴───────────┴───────────┘
Get the covariance matrix
matrix(m)Get the correlation matrix
matrix(m, correlation=true)Call asymmetric uncertainty intervals ("Minos errors")
Call minos! to get parameter errors
minos!(m)┌──────────────┬──────────────┬───────────┬────────────┬──────────────┐
│ FCN │ Method │ Ncalls │ Iterations │ Elapsed time │
│ 9.033 │ migrad │ 119 │ 8 │ 481.346 µs │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Valid Min. │ Valid Param. │ Above EDM │ Call limit │ Edm │
│ true │ true │ false │ false │ 7.48753e-19 │
├──────────────┼──────────────┼───────────┼────────────┼──────────────┤
│ Hesse failed │ Has cov. │ Accurate │ Pos. def. │ Forced │
│ false │ true │ true │ true │ false │
└──────────────┴──────────────┴───────────┴────────────┴──────────────┘
┌───┬──────┬──────────┬─────────────┬────────────┬───────────┬────────┬────────┬───────┬───────┐
│ │ Name │ Value │ Hesse Error │ Minos- │ Minos+ │ Limit- │ Limit+ │ Fixed │ Const │
├───┼──────┼──────────┼─────────────┼────────────┼───────────┼────────┼────────┼───────┼───────┤
│ 1 │ a │ 0.911937 │ 0.0907645 │ -0.0902591 │ 0.0902591 │ │ │ │ │
│ 2 │ b │ 2.73544 │ 0.831558 │ -0.82689 │ 0.82689 │ │ │ │ │
│ 3 │ c │ -1.50305 │ 1.99903 │ -1.98779 │ 1.98779 │ │ │ │ │
│ 4 │ d │ 0.765475 │ 1.31179 │ -1.30442 │ 1.30442 │ │ │ │ │
└───┴──────┴──────────┴─────────────┴────────────┴───────────┴────────┴────────┴───────┴───────┘
┌───┬───────────┬───────────┬───────────┬───────────┐
│ │ a │ b │ c │ d │
├───┼───────────┼───────────┼───────────┼───────────┤
│ a │ 1.0 │ -0.757518 │ 0.593245 │ -0.499613 │
│ b │ -0.757518 │ 1.0 │ -0.958513 │ 0.89783 │
│ c │ 0.593245 │ -0.958513 │ 1.0 │ -0.984329 │
│ d │ -0.499613 │ 0.89783 │ -0.984329 │ 1.0 │
└───┴───────────┴───────────┴───────────┴───────────┘
Get the asymmetric errors
m.minos |> showUnion{Nothing, Minuit2.ROOT!Minuit2!MinosError}[┌──────────┬────────────┬───────────┐
│ 1 │ valid │ │
├──────────┼────────────┼───────────┤
│ Error │ -0.0902591 │ 0.0902591 │
│ Valid │ true │ true │
│ At Limit │ false │ false │
│ Max Fcn │ false │ false │
│ New Min │ false │ false │
└──────────┴────────────┴───────────┘
, ┌──────────┬──────────┬─────────┐
│ 2 │ valid │ │
├──────────┼──────────┼─────────┤
│ Error │ -0.82689 │ 0.82689 │
│ Valid │ true │ true │
│ At Limit │ false │ false │
│ Max Fcn │ false │ false │
│ New Min │ false │ false │
└──────────┴──────────┴─────────┘
, ┌──────────┬──────────┬─────────┐
│ 3 │ valid │ │
├──────────┼──────────┼─────────┤
│ Error │ -1.98779 │ 1.98779 │
│ Valid │ true │ true │
│ At Limit │ false │ false │
│ Max Fcn │ false │ false │
│ New Min │ false │ false │
└──────────┴──────────┴─────────┘
, ┌──────────┬──────────┬─────────┐
│ 4 │ valid │ │
├──────────┼──────────┼─────────┤
│ Error │ -1.30442 │ 1.30442 │
│ Valid │ true │ true │
│ At Limit │ false │ false │
│ Max Fcn │ false │ false │
│ New Min │ false │ false │
└──────────┴──────────┴─────────┘
]Minos can fail, check messages:
- "Valid": everything is chipper
- "At Limit": Minos hit parameter limit before finishing contour
- "Max FCN": Minos reached call limit before finishing contour
- "New Min": Minos found a new minimum while scanning
get the minos errors
d = m.minos["d"]
d.lower, d.upper(-1.3044247452958948, 1.3044247449130866)Plot the parameters with the asymmetric errors plot parameters with errors
v = m.values |> collect
ve = m.errors |> collect
vm = m.minos |> values
vmean = [(e.upper - e.lower) / 2 for e in vm]
npar = length(v)
indices = 1:npar
plot(indices .- 0.05, v, yerr=ve, seriestype=:scatter, label="Hesse errors")
plot!(indices .+ 0.05, v, yerr=vmean, seriestype=:scatter, label="Minos errors")Built-in plotting
These are some functions that provide some plotting functionality. There are implemented as an extension module to Minuit2, which is loaded when the module Plots is available.
draw_contour(m, "d", "c", bound=4, subtract_min=false)Draw the contour plot for the parameters d and c
draw_mncontour(m, "d", "c", cl=1:4)Draw the profile plot for the parameter d
draw_profile(m, "d", bound=2)Draw the profile plot for the parameter d
draw_mnprofile(m, "d", size=20)Summary
In this example we have shown how to use the Minuit2 module to fit a line model to some data and how to estimate the parameter uncertainties. We have also shown how to fit a polynomial model with a flexible number of parameters.
This page was generated using Literate.jl.