probfit Basic Tutorial¶

#lets make s straight line with gaussian(mu=0,sigma=1) noise
npr.seed(0)
x = linspace(0,10,20) 
y = 3*x+15+ randn(20)
err = np.array([1]*20)
errorbar(x,y,err,fmt='.');

describe(line)

#Chi^2 regression is just a callable object nothing special about it
describe(chi2)

#minimize it
#yes it gives you a heads up that you didn't give it initial value
#we can ignore it for now
minimizer = iminuit.Minuit(chi2) #see iminuit tutorial on how to give initial value/range/error
minimizer.migrad(); #very stable robust minimizer
#you can look at your terminal to see what it is doing;

#lets see our results
print minimizer.values
print minimizer.errors

#or a pretty printing
minimizer.print_fmin()

#let's visualize our line
chi2.draw(minimizer)
#looks good

#Sometime we want the error matrix
print 'error matrix:\n', minimizer.matrix()
#or correlation matrix
print 'correlation matrix:\n', minimizer.matrix(correlation=True)
#or a pretty html
#despite the method named error_matrix it's actually a correlation matrix
minimizer.print_matrix()

#First lets make our data
npr.seed(0)
data = randn(10000)*4+1
#sigma = 4 and mean = 1
hist(data,bins=100,histtype='step');

#Let's fit
minimizer = iminuit.Minuit(cost,sigma=3.) #notice here that we give initial value to sigma
#but most of the time we want to see it before fitting
cost.draw(minimizer)

minimizer.migrad() #very stable minimization algorithm
#like in all binned fit with long zero tail. It will have to do something about the zero bin
#dist_fit handle them gracefully but will give you a head up

#you can see if your fit make any sense too
cost.draw(minimizer)#uncertainty is given by symetric poisson

#let's see the result
print 'Value:', minimizer.values
print 'Error:', minimizer.errors

#That printout can get out of hand quickly
minimizer.print_fmin()
#and correlation matrix
#will not display well in firefox(tell them to fix writing-mode:)
minimizer.print_matrix()

#how about making sure the error making sense
minimizer.draw_mnprofile('mu');

#2d contour error
#you can notice that it takes sometime to draw
#we will this is because our PDF is defined in Python
#we will show how to speed this up later
minimizer.draw_mncontour('mu','sigma');

#we will use the same data as in the previous example
npr.seed(0)
data = npr.randn(10000)*4+1
#sigma = 4 and mean = 1
hist(data, bins=100, histtype='step');

#very useful method to look at function signature
describe(extended_pdf) #you can look at what your pdf means

#Chi^2 distribution fit is really bad for distribution with long tail
#since when bin count=0... poisson error=0 and blows up chi^2
#so give it some range
chi2 = BinnedChi2(extended_pdf,data,bound=(-7,10))
minimizer = iminuit.Minuit(chi2,sigma=1)
minimizer.migrad()

chi2.draw(minimizer)

#and usual deal
minimizer.print_fmin()
minimizer.print_matrix() 

#same data
npr.seed(0)
data = npr.randn(10000)*4+1
#sigma = 4 and mean = 1
hist(data,bins=100,histtype='step');

#cost function 
ublh = UnbinnedLH(cython_PDF,data)
minimizer = iminuit.Minuit(ublh,sigma=2.)
minimizer.set_up(0.5)#remember this is likelihood
minimizer.migrad()#yes amazingly fast
ublh.show(minimizer)
minimizer.print_fmin()
minimizer.print_matrix() 

#remember how slow it was?
#now it's super fast(and it's even unbinned likelihood)
minimizer.draw_mnprofile('mu');

import probfit.pdf
print dir(probfit.pdf)
print describe(gaussian)
print type(gaussian)

ublh = UnbinnedLH(gaussian,data)
minimizer = iminuit.Minuit(ublh,sigma=2.)
minimizer.set_up(0.5)#remember this is likelihood
minimizer.migrad()#yes amazingly fast
ublh.draw(minimizer, show_errbars='normal') # control how fit is displayed too

#lets first generate a crystal ball sample
#dist_fit has builtin toy generation capability
#lets introduce crystal ball function
#http://en.wikipedia.org/wiki/Crystal_Ball_function
#it's simply gaussian with power law tail
#normally found in energy deposited in crystals
#impossible to normalize analytically
#and normalization will depend on shape parameters
describe(crystalball)

npr.seed(0)
bound = (-1,2)
data = gen_toy(crystalball,10000,bound=bound,alpha=1.,n=2.,mean=1.,sigma=0.3,quiet=False)
#quiet = False tells it to plot out original function
#toy histogram and poisson error from both orignal distribution and toy

#To fit this we need to tell normalized our crystal ball PDF
#this is done with trapezoid rule with simple cache mechanism
#can be done by Normalized functor
ncball = Normalized(crystalball,bound)
#this can also bedone with declerator
#@normalized_function(bound)
#def myPDF(x,blah):
#    return blah
print 'unnorm:', crystalball(1.0,1,2,1,0.3)
print '  norm:', ncball(1.0,1,2,1,0.3)

#it has the same signature as the crystalball
describe(ncball)

#now we can fit as usual
ublh = UnbinnedLH(ncball,data)
minimizer = iminuit.Minuit(ublh,
    alpha=1.,n=2.1,mean=1.2,sigma=0.3)
minimizer.set_up(0.5)#remember this is likelihood
minimizer.migrad()#yes amazingly fast Normalize is written in cython
ublh.show(minimizer)
#crystalball function is nortorious for its sensitivity on n parameter
#dist_fit give you a heads up where it might have float overflow

from probfit import integrate1d
def line(x, m, c):
    return m*x+c

#compute integral of line from x=(0,1) using 10 intevals with m=1. and c=2.
#all probfit internal use this
print integrate1d(line, (0,1), 10, (1.,2.) ) # no integrate method available probfit use simpson3/8

#Let us illustrate the point by forcing it to have integral that's off by
#factor of two
def wrong_line_integral(bound, nint, m, c):
    a, b = bound
    return 2*(m*(b**2/2.-a**2/2.)+c*(b-a)) # I know this is wrong 
line.integrate = wrong_line_integral
# line.integrate = lambda bound, nint, m, c: blah blah # this works too
print integrate1d(line, (0,1), 10, (1.,2.))

#crystalball is nortoriously sensitive to initial parameter
#now it is a good time to show what happen when things...doesn't fit
ublh = UnbinnedLH(ncball,data)
minimizer = iminuit.Minuit(ublh)#NO initial value
minimizer.set_up(0.5)#remember this is likelihood
minimizer.migrad()#yes amazingly fast but tons of read
#Remember there is a heads up

ublh.show(minimizer)#it bounds to fails

minimizer.migrad_ok(), minimizer.matrix_accurate()
#checking these two method give you a good sign

#fix it by giving it initial value/error/limit or fixing parameter see iminuit Tutorial
#now we can fit as usual
ublh = UnbinnedLH(ncball,data)
minimizer = iminuit.Minuit(ublh,
    alpha=1.,n=2.1,mean=1.2,sigma=0.3)
minimizer.set_up(0.5)#remember this is likelihood
minimizer.migrad()#yes amazingly fast. Normalize is written in cython
ublh.show(minimizer)

besttry = try_uml(ncball,data,alpha=1.,n=2.1,mean=1.2,sigma=0.3)

#or you can try multiple
#too many will just confuse you
besttry = try_uml(ncball,data,alpha=1.,n=2.1,mean=[1.2,1.1],sigma=[0.3,0.5])
print besttry #and you can find which one give you minimal unbinned likelihood

peak1 = npr.randn(3000)*0.2
peak2 = npr.randn(5000)*0.1+4
bg = gen_toy(lambda x : (x+2)**2, 20000,(-2,5))
all_data = np.concatenate([peak1,peak2,bg])
hist((peak1,peak2,bg,all_data),bins=200,histtype='step',range=(-2,5));

%load_ext cythonmagic

%%cython
cimport cython
from probfit import Normalized, gaussian

@cython.binding(True)
def poly(double x,double a,double b, double c):
    return a*x*x+b*x+c

#remember linear function is not normalizable for -inf ... +inf
nlin = Normalized(poly,(-1,5))

#our extended PDF for 3 types of signal
@cython.binding(True)
def myPDF(double x, 
    double a, double b, double c, double nbg,
    double mu1, double sigma1, double nsig1,
    double mu2, double sigma2, double nsig2):

    cdef double NBG = nbg*nlin(x,a,b,c)
    cdef double NSIG1 = nsig1*gaussian(x,mu1,sigma1) 
    cdef double NSIG2 = nsig2*gaussian(x,mu2,sigma2)
    return NBG + NSIG1 + NSIG2

print describe(myPDF)

#lets use what we just learned
#for complicated function good initial value(and initial step) is crucial
#if it doesn't converge try play with initial value and initial stepsize(error_xxx)
besttry = try_binlh(myPDF,all_data, 
    a=1.,b=2.,c=4.,nbg=20000.,
    mu1=0.1,sigma1=0.2,nsig1=3000.,
    mu2=3.9,sigma2=0.1,nsig2=5000.,extended=True, bins=300, bound=(-1,5) )
print besttry

binlh = BinnedLH(myPDF,all_data, bins=200, extended=True, bound=(-1,5))
#too lazy to type initial values from what we try?
#use keyword expansion **besttry
#need to put in initial step size for mu1 and mu2 
#with error_mu* otherwise it won't converge(try it yourself)
minimizer = iminuit.Minuit(binlh, error_mu1=0.1, error_mu2=0.1, **besttry)

minimizer.migrad();

binlh.show(minimizer)

minimizer.print_fmin()

npr.seed(0)
#Lets make two gaussian
data1 = npr.randn(10000)+3
data2 = npr.randn(10000)-2
hist([data1,data2],bins=100,histtype='step',label=['data1','data2']);

#remember this is nothing special about builtin cost function
#except some utility function like draw and show
ulh1 = UnbinnedLH(gaussian,data1)
ulh2 = UnbinnedLH(gaussian,data2)
print describe(ulh1)
print describe(ulh2)

minimizer = iminuit.Minuit(simul_lh,sigma=0.5)
minimizer.set_up(0.5)#remember it's likelihood
minimizer.migrad();

minimizer.print_fmin()
minimizer.print_matrix()
results = minimizer.values

draw_compare_hist(gaussian,[results['mu1'],results['sigma']],data1,normed=True);
draw_compare_hist(gaussian,[results['mu2'],results['sigma']],data2,normed=True);