In [1]:
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

Really Quick Start

Let go through a quick course about how to minimize things. If you use PyMinuit before you will find that iminuit is very similar to PyMinuit.

Our first example is about trying to minimize a simple function:

$f(x,y,z) = (x-1)^2 + (y-2)^2 + (z-3)^2 - 1$

We know easily that the answer has is $x=1$, $y=2$, $z=3$ and the minimum value should be $-1$

In [2]:
def f(x,y,z):
    return (x-1.)**2 + (y-2*x)**2 + (z-3.*x)**2 -1.

iminuit relies on Python introspection. If you wonder what kind of function iminuit sees, you can use

In [3]:
describe(f) 
Out[3]:
['x', 'y', 'z']

Construct Minuit Object

To minimize we need to construct a Minuit object

In [4]:
m=Minuit(f, x=2, error_x=0.2, limit_x=(-10.,10.), y=3., fix_y=True, print_level=1)
-c:1: InitialParamWarning: errordef is not given. Default to 1.
-c:1: InitialParamWarning: Parameter z does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter z is floating but does not have initial step size. Assume 1.

The initial value/error are optional but it's nice to do it and here is how to use it

  • x=2 set intial value of x to 2
  • error_x=0.2 set the initial stepsize
  • limit_x = (-1,1) set the range for x
  • y=2, fix_y=True fix y value to 2

We did not put any constain on z. Minuit will howerver warn you about missig initial error/step(using python builtin warning).

Check that initial setting you put in is what you want

In [5]:
m.print_param()#or call print_initial_param
#bonus: if you click the + button on the top left corner it will show latex table
#which you can copy and paste to your beamer/latex document
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 2.000000e+00 2.000000e-01 0.000000e+00 0.000000e+00 -10.0 10.0
2 y 3.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 FIXED
3 z 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00

Run Migrad

Migrad performs Variable Metric Minimization. In a nutshell, it combines steepest descends algorithm along with line search strategy. Migrad is very popular in high energy physics field because of its robustness.

In [6]:
#Minimize
m.migrad();
#notice also in your prompt it prints out progress

FCN = -0.799999984093 NFCN = 36 NCALLS = 36
EDM = 1.5906866957e-08 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 1.399979e+00 4.470607e-01 0.000000e+00 0.000000e+00 -10.0 10.0
2 y 3.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 FIXED
3 z 4.199819e+00 1.673317e+00 0.000000e+00 0.000000e+00

Out[6]:
({'hesse_failed': False, 'has_reached_call_limit': False, 'has_accurate_covar': True, 'has_posdef_covar': True, 'up': 1.0, 'edm': 1.590686695699253e-08, 'is_valid': True, 'is_above_max_edm': False, 'has_covariance': True, 'has_made_posdef_covar': False, 'has_valid_parameters': True, 'fval': -0.7999999840931111, 'nfcn': 36},
 [{'is_const': False, 'name': 'x', 'has_limits': True, 'value': 1.3999787801966246, 'number': 0L, 'has_lower_limit': True, 'upper_limit': 10.0, 'lower_limit': -10.0, 'has_upper_limit': True, 'error': 0.447060710371745, 'is_fixed': False},
  {'is_const': False, 'name': 'y', 'has_limits': False, 'value': 3.0, 'number': 1L, 'has_lower_limit': False, 'upper_limit': 0.0, 'lower_limit': 0.0, 'has_upper_limit': False, 'error': 1.0, 'is_fixed': True},
  {'is_const': False, 'name': 'z', 'has_limits': False, 'value': 4.199819483888685, 'number': 2L, 'has_lower_limit': False, 'upper_limit': 0.0, 'lower_limit': 0.0, 'has_upper_limit': False, 'error': 1.6733166635386176, 'is_fixed': False}])

migrad summary table give you a nice overview of fit status.

  • All blocks should be green.
  • Red means something bad.
  • Yellow is a caution(fit is good but you ran over call limit)

You can use the return value of migrad() to check fit status. Most important field is is_valid.

Accessing values and errors

Accessing Values

In [7]:
#and this is how you get the the value
print('parameters', m.parameters)
print('args', m.args)
print('value', m.values)
parameters ['x', 'y', 'z']
args (1.3999787801966246, 3.0, 4.199819483888685)
value {'y': 3.0, 'x': 1.3999787801966246, 'z': 4.199819483888685}

Error(parabolic)

In [8]:
#and the error
print('error', m.errors)
error {'y': 1.0, 'x': 0.447060710371745, 'z': 1.6733166635386176}

Function minimum

In [9]:
#and function value at the minimum
print('fval', m.fval)
#Tip: you can also obtain value at current state by
print('current state', f(*m.args))
fval -0.799999984093
current state -0.799999984093

Correlation and Covariance Matrix

In [10]:
#covariance, correlation matrix
#remember y is fixed
print('covariance', m.covariance)
print('matrix()', m.matrix()) #covariance
print('matrix(correlation=True)', m.matrix(correlation=True)) #correlation
m.print_matrix() #correlation
#again click the + button on the top left corner for latex code
covariance {('z', 'z'): 2.7999886564760113, ('x', 'z'): 0.5999969687373324, ('x', 'x'): 0.1999992395679503, ('z', 'x'): 0.5999969687373324}
matrix() ((0.1999992395679503, 0.5999969687373324), (0.5999969687373324, 2.7999886564760113))
matrix(correlation=True) ((1.0, 0.8017828234103097), (0.8017828234103097, 1.0))
+
x
z
x 1.00 0.80
z 0.80 1.00

Fit status

In [11]:
#get mimization status
print(m.get_fmin())
print(m.get_fmin().is_valid)
{'hesse_failed': False, 'has_reached_call_limit': False, 'has_accurate_covar': True, 'has_posdef_covar': True, 'up': 1.0, 'edm': 1.590686695699253e-08, 'is_valid': True, 'is_above_max_edm': False, 'has_covariance': True, 'has_made_posdef_covar': False, 'has_valid_parameters': True, 'fval': -0.7999999840931111, 'nfcn': 36}
True

Contour and $\chi^2$/Likelihood profile

$\chi^2$ and contour can be obtained easily

In [12]:
#minos contour
xminos, zminos, ctr = m.mncontour('x','z')
fill(*zip(*ctr)) #looks kinda ugly, right?
Out[12]:
[<matplotlib.patches.Polygon at 0x104026f10>]
In [13]:
#drawing it nicely
m.draw_mncontour('x','z', nsigma=4);
In [14]:
#or you can get the gridded data
x, y, g, r = m.mncontour_grid('x','z', nsigma=3) # r is the raw data
pcolormesh(x,y,g)
colorbar()
Out[14]:
<matplotlib.colorbar.Colorbar instance at 0x1042d4050>
In [15]:
#1D value Scan
x,y = m.profile('x',subtract_min=True);
plot(x,y) #if you have matplotlib
Out[15]:
[<matplotlib.lines.Line2D at 0x10485b1d0>]
In [16]:
#we also provide convenience wrapper for drawing it
m.draw_profile('x');
In [17]:
#2d contour NOT minos contour
x,y,z = m.contour('x','y',subtract_min=True)
cs = contour(x,y,z)
clabel(cs)
Out[17]:
<a list of 8 text.Text objects>
In [18]:
#or a convenience wrapper
m.draw_contour('x','z');

Hesse and Minos

Hesse

Hesse find the error by finding the inverse of second derivative matrix(hessian). The error assume parabolic shape at the minimum. Hesse error is symmetric by construct. Hesse is always called at the end of migrad to get the error. You normally don't need to call it manually.

In [19]:
m.hesse()
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 1.399979e+00 4.470591e-01 0.000000e+00 0.000000e+00 -10.0 10.0
2 y 3.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 FIXED
3 z 4.199819e+00 1.673311e+00 0.000000e+00 0.000000e+00
+
x
z
x 1.00 0.80
z 0.80 1.00
Out[19]:
[{'is_const': False, 'name': 'x', 'has_limits': True, 'value': 1.3999787801966246, 'number': 0L, 'has_lower_limit': True, 'upper_limit': 10.0, 'lower_limit': -10.0, 'has_upper_limit': True, 'error': 0.4470591480557786, 'is_fixed': False},
 {'is_const': False, 'name': 'y', 'has_limits': False, 'value': 3.0, 'number': 1L, 'has_lower_limit': False, 'upper_limit': 0.0, 'lower_limit': 0.0, 'has_upper_limit': False, 'error': 1.0, 'is_fixed': True},
 {'is_const': False, 'name': 'z', 'has_limits': False, 'value': 4.199819483888685, 'number': 2L, 'has_lower_limit': False, 'upper_limit': 0.0, 'lower_limit': 0.0, 'has_upper_limit': False, 'error': 1.673310811142077, 'is_fixed': False}]

minos

minos multidimensionally scan likelihood/$\chi^2$ until to find the contour where the value of the cost function increase by UP(see set_up). It takes really long time but give the correct error(unless it fails).

In [20]:
m.minos() #call m.minos('x') if you need minos error for just 1 variable
print(m.get_merrors()['x'])
print(m.get_merrors()['x'].lower)
print(m.get_merrors()['x'].upper)
Minos status for x: VALID
Error -0.447192373006 0.447234827127
Valid True True
At Limit False False
Max FCN False False
New Min False False
Minos status for z: VALID
Error -1.67313956898 1.67350060134
Valid True True
At Limit False False
Max FCN False False
New Min False False
{'lower_new_min': False, 'upper': 0.44723482712747076, 'lower': -0.4471923730061797, 'at_lower_limit': False, 'min': 1.3999787801966246, 'at_lower_max_fcn': False, 'is_valid': True, 'upper_new_min': False, 'at_upper_limit': False, 'lower_valid': True, 'upper_valid': True, 'at_upper_max_fcn': False, 'nfcn': 16L}
-0.447192373006
0.447234827127

Printing Out Nice Tables

you can force use print_* to do various html display

In [21]:
m.print_param()
m.print_matrix()
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 1.399979e+00 4.470591e-01 -4.471924e-01 4.472348e-01 -10.0 10.0
2 y 3.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00 FIXED
3 z 4.199819e+00 1.673311e+00 -1.673140e+00 1.673501e+00
+
x
z
x 1.00 0.80
z 0.80 1.00

Alternative Ways to define function

Cython

If you want speed with minimal code change this is the way to do it. This is a quick way to use cython. For a hard core cython see hard-core-tutorial.ipynb.

In [22]:
#sometimes we want speeeeeeed
%load_ext Cython
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [23]:
%%cython --force
cimport cython
import numpy as np
cimport numpy as np #overwritten those from python with cython

[email protected](True) this works too
@cython.embedsignature(True)#dump signature in pydoc so describe can extract signature
def cython_f(double x,double y,double z):
    return (x-1.)**2 + (y-2.)**2 + (z-3.)**2 -1.
In [24]:
#you can always see what iminuit will see
print(describe(cython_f))
['x', 'y', 'z']
In [25]:
m = Minuit(cython_f)
m.migrad()
print(m.values)
-c:1: InitialParamWarning: errordef is not given. Default to 1.
-c:1: InitialParamWarning: Parameter x does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter x is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter y does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter y is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter z does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter z is floating but does not have initial step size. Assume 1.

FCN = -1.0 NFCN = 36 NCALLS = 36
EDM = 8.24353485149e-23 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 1.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00
2 y 2.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00
3 z 3.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00

{'y': 2.0000000000047877, 'x': 1.000000000002637, 'z': 3.000000000007182}

Callable object ie: call

This is useful if you need to bound your object to some data

In [26]:
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [27]:
x = [1,2,3,4,5]
y = [2,4,6,8,10]# y=2x
class StraightLineChi2:
    def __init__(self,x,y):
        self.x = x
        self.y = y
    def __call__(self,m,c): #lets try to find slope and intercept
        chi2 = sum((y - m*x+c)**2 for x,y in zip(self.x,self.y))
        return chi2
In [28]:
chi2 = StraightLineChi2(x,y)
describe(chi2)
Out[28]:
['m', 'c']
In [29]:
m = Minuit(chi2)
m.migrad()
print(m.values)
-c:1: InitialParamWarning: errordef is not given. Default to 1.
-c:1: InitialParamWarning: Parameter m does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter m is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter c does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter c is floating but does not have initial step size. Assume 1.

FCN = 6.31477238173e-26 NFCN = 32 NCALLS = 32
EDM = 6.32720690373e-26 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 m 2.000000e+00 3.162278e-01 0.000000e+00 0.000000e+00
2 c -4.884981e-14 1.048809e+00 0.000000e+00 0.000000e+00

{'c': -4.884981308350689e-14, 'm': 2.00000000000002}

Faking a function signature

This is missing from PyMinuit. iminuit allows you to take funciton sinature by using func_code.co_varnames and func_code.co_argcount. This is very useful for making a higher order function that takes PDF and data in to calculate appropriate cost function.

In [30]:
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [31]:
#this is very useful if you want to build a generic cost functor
#this is actually how probfit is implemented
from iminuit.util import make_func_code
x = [1,2,3,4,5]
y = [2,4,6,8,10]# y=2x
class Chi2Functor:
    def __init__(self,f,x,y):
        self.f = f
        self.x = x
        self.y = y
        f_sig = describe(f)
        #this is how you fake function 
        #signature dynamically
        self.func_code = make_func_code(f_sig[1:])#docking off independent variable
        self.func_defaults = None #this keeps np.vectorize happy
    def __call__(self,*arg):
        #notice that it accept variable length
        #positional arguments
        chi2 = sum((y-self.f(x,*arg))**2 for x,y in zip(self.x,self.y))
        return chi2
In [32]:
def linear(x,m,c):
    return m*x+c

def parabola(x,a,b,c):
    return a*x**2 + b*x + c 
In [33]:
linear_chi2 = Chi2Functor(linear,x,y)
describe(linear_chi2)
Out[33]:
['m', 'c']
In [34]:
m = Minuit(linear_chi2)
m.migrad();
print(m.values)
-c:1: InitialParamWarning: errordef is not given. Default to 1.
-c:1: InitialParamWarning: Parameter m does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter m is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter c does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter c is floating but does not have initial step size. Assume 1.

FCN = 6.48916980635e-26 NFCN = 32 NCALLS = 32
EDM = 6.48992044656e-26 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 m 2.000000e+00 3.162278e-01 0.000000e+00 0.000000e+00
2 c 3.375078e-14 1.048809e+00 0.000000e+00 0.000000e+00

{'c': 3.375077994860476e-14, 'm': 2.000000000000025}
In [35]:
#now here is the beauty
#you can use the same Chi2Functor now for parabola
parab_chi2 = Chi2Functor(parabola,x,y)
describe(parab_chi2)
Out[35]:
['a', 'b', 'c']
In [36]:
m = Minuit(parab_chi2,x,y)
m.migrad()
print(m.values)
-c:1: InitialParamWarning: errordef is not given. Default to 1.
-c:1: InitialParamWarning: Parameter a does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter a is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter b does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter b is floating but does not have initial step size. Assume 1.
-c:1: InitialParamWarning: Parameter c does not have initial value. Assume 0.
-c:1: InitialParamWarning: Parameter c is floating but does not have initial step size. Assume 1.

FCN = 1.50034031957e-18 NFCN = 59 NCALLS = 59
EDM = 1.50033920511e-18 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 a -4.761125e-11 2.672612e-01 0.000000e+00 0.000000e+00
2 b 2.000000e+00 1.634451e+00 0.000000e+00 0.000000e+00
3 c 8.948271e-10 2.144761e+00 0.000000e+00 0.000000e+00

{'a': -4.761124827723506e-11, 'c': 8.948271013053954e-10, 'b': 1.9999999999086318}

Last Resort: Forcing function signature

built-in function normally do not have signature. Function from swig belongs in this categories. Python intro spection will fails and we have to force function signature.

In [37]:
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [38]:
%%cython
#sometimes you get a function with absolutely no signature
def nosig_f(x,y):
    return x**2+(y-4)**2
In [39]:
#something from swig will give you a function with no
#signature information
try:
    describe(nosig_f)#it raise error
except Exception as e:
    print(e)
Unable to obtain function signature
In [40]:
#Use forced_parameters
m = Minuit(nosig_f, forced_parameters=('x','y'))
-c:2: InitialParamWarning: errordef is not given. Default to 1.
-c:2: InitialParamWarning: Parameter x does not have initial value. Assume 0.
-c:2: InitialParamWarning: Parameter x is floating but does not have initial step size. Assume 1.
-c:2: InitialParamWarning: Parameter y does not have initial value. Assume 0.
-c:2: InitialParamWarning: Parameter y is floating but does not have initial step size. Assume 1.
In [41]:
m.migrad()
print(m.values)

FCN = 9.31236296332e-23 NFCN = 24 NCALLS = 24
EDM = 9.31236296802e-23 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Parab Error Minos Error- Minos Error+ Limit- Limit+ FIXED
1 x 0.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00
2 y 4.000000e+00 1.000000e+00 0.000000e+00 0.000000e+00

{'y': 4.00000000000965, 'x': 0.0}
In [42]:
print(f)
<built-in method f of mtrand.RandomState object at 0x1002ab4e0>

Frontend

Frontend affects how the output from migrad/minos etc are displayed. iminuit is shipped with two frontends. ConsoleFrontend print in text format and HtmlFrontend print html object to Ipython notebook. When you construct Minuit object the front end is selected automatically. If you are in IPython it will use Html frontend; otherwise, it will use console fronend. You can force Minuit to use frontend of your choice too.

In [43]:
%pylab inline
from iminuit import Minuit, describe, Struct
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
In [44]:
#this is just showing off console frontend (you can force it)
from iminuit.frontends import ConsoleFrontend
def f(x,y,z):
    return (x-1.)**2 + (y-2.)**2 + (z-3.)**2 -1.
m = Minuit(f, frontend=ConsoleFrontend())
-c:5: InitialParamWarning: errordef is not given. Default to 1.
-c:5: InitialParamWarning: Parameter x does not have initial value. Assume 0.
-c:5: InitialParamWarning: Parameter x is floating but does not have initial step size. Assume 1.
-c:5: InitialParamWarning: Parameter y does not have initial value. Assume 0.
-c:5: InitialParamWarning: Parameter y is floating but does not have initial step size. Assume 1.
-c:5: InitialParamWarning: Parameter z does not have initial value. Assume 0.
-c:5: InitialParamWarning: Parameter z is floating but does not have initial step size. Assume 1.
In [45]:
m.migrad();
**************************************************
*                     MIGRAD                     *
**************************************************

**********************************************************************
---------------------------------------------------------------------------------------
fval = -1.0 | nfcn = 36 | ncalls = 36
edm = 8.243534851490919e-23 (Goal: 1e-05) | up = 1.0
---------------------------------------------------------------------------------------
|          Valid |    Valid Param | Accurate Covar |         Posdef |    Made Posdef |
---------------------------------------------------------------------------------------
|           True |           True |           True |           True |          False |
---------------------------------------------------------------------------------------
|     Hesse Fail |        Has Cov |      Above EDM |                |  Reach calllim |
---------------------------------------------------------------------------------------
|          False |           True |          False |             '' |          False |
---------------------------------------------------------------------------------------

----------------------------------------------------------------------------------------------
|      | Name  |  Value   | Para Err |   Err-   |   Err+   |  Limit-  |  Limit+  |          |
----------------------------------------------------------------------------------------------
|    0 |     x =  1       |  1       |          |          |          |          |          |
|    1 |     y =  2       |  1       |          |          |          |          |          |
|    2 |     z =  3       |  1       |          |          |          |          |          |
----------------------------------------------------------------------------------------------

**********************************************************************
In [ ]: