Using external minimizers

We show how to use an external minimizer to find the minimum of a function and then use iminuit to compute the parameter uncertainties.

We will demonstrate this with a maximum-likelihood fit of a normal distribution, which is carried out with scipy.optimize.minimize. iminuit is then used to compute the parameter uncertainties.

Note: iminuit can call the scipy minimizers directly with Minuit.scipy, so scipy.optimize.minimize is only used here to demonstrate the general approach.

from iminuit import Minuit
import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize

# normally distributed data
rng = np.random.default_rng(1)
x = rng.normal(size=1000)


# negative log-likelihood for a normal distribution
def nll(par):
    return -np.sum(norm.logpdf(x, par[0], par[1]))


nll.errordef = Minuit.LIKELIHOOD

# minimize nll with scipy.optimize.minimize
result = minimize(nll, np.ones(2))
result

  message: Desired error not necessarily achieved due to precision loss.
  success: False
   status: 2
      fun: 1405.1044194490182
        x: [-5.425e-02  9.863e-01]
      nit: 12
      jac: [ 1.526e-05  0.000e+00]
 hess_inv: [[ 5.917e-04  7.100e-05]
            [ 7.100e-05  4.437e-04]]
     nfev: 48
     njev: 16

# initialize Minuit with the fit result from scipy.optimize.minimize
m = Minuit(nll, result.x)
m.hesse()  # this also works without calling MIGRAD before

External
FCN = 1405	Nfcn = 25
EDM = 9.78e-14 (Goal: 0.0001)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	-0.054	0.031
1	x1	0.986	0.022

	x0	x1
x0	0.000973	-0
x1	-0	0.000486

We can also compute the “Hesse errors” at any other point than the minimum. These cannot be interpreted as parameter uncertainties, they are just some numbers related to the second derivative of the cost function at that point.

m.values = (1.0, 0.5)
m.hesse()

External
FCN = 4394	Nfcn = 48
EDM = 2.42e+03 (Goal: 0.0001)
INVALID Minimum	ABOVE EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	1.000	0.031
1	x1	0.500	0.006

	x0	x1
x0	0.000964	0.17e-3 (0.861)
x1	0.17e-3 (0.861)	4.01e-05

Minuit now reports that the minimum is invalid, which is correct, but it does not matter for the Hesse errors, which are computed anyway.

Likewise, it one can also run MINOS to get MINOS estimates. Note that MINOS can fail if the starting point is not actually a minimum. So here we reset the values to the solution found by scipy.optimize.

m.values = result.x
m.minos()

External
FCN = 1405	Nfcn = 111
EDM = 9.78e-14 (Goal: 0.0001)
Valid Minimum	Below EDM threshold (goal x 10)
No parameters at limit	Below call limit
Hesse ok	Covariance accurate

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	x0	-0.054	0.031	-0.031	0.031
1	x1	0.986	0.022	-0.022	0.022

	x0		x1
Error	-0.031	0.031	-0.022	0.022
Valid	True	True	True	True
At Limit	False	False	False	False
Max FCN	False	False	False	False
New Min	False	False	False	False

	x0	x1
x0	0.000973	-0
x1	-0	0.000486

We can see that MINOS ran successfully. The Hesse Errors were also updated, because MINOS needs HESSE to run first. HESSE is called automatically in this case.