How to do model fitting in Python

Curve fitting - How?¶

For the purpose of illustrating an example, I will introduce a limb model which I have been working on during my PhD.¶

In [1]:
## importing the libraries
import numpy as np
import math as m
import scipy
import matplotlib.pyplot as plt
from scipy import special

In [2]:
def Erfc(x,sigma):

y = special.erfc(x/(sigma*np.sqrt(2)))
return y

## Introducing the model to be used later for the fitting

def SL_fit(x,w1,w2,w3,s1,s2,s3):

f = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+ (1-w1-w2-w3))
return f


The first method we will use is via scipy.optimize.curve_fit¶

In [4]:
from scipy.optimize import curve_fit

In [5]:
#This is the jacobian matrx, which is the derivative of the fitting function with respect to every free parameter
def DF(x,w1,w2,w3,s1,s2,s3):
dfw1 = 0.5*(Erfc(x,s1)-1)
dfw2 = 0.5*(Erfc(x,s2)-1)
dfw3 =  0.5*(Erfc(x,s3)-1)
dfs1 = w1*x*np.exp(-x**2/(2*s1**2))/ (np.sqrt(2*np.pi)*s1**2)
dfs2 = w2*x*np.exp(-x**2/(2*s2**2))/ (np.sqrt(2*np.pi)*s2**2)
dfs3 = w3*x*np.exp(-x**2/(2*s3**2))/ (np.sqrt(2*np.pi)*s3**2)
return np.transpose(np.array([dfw1,dfw2,dfw3,dfs1,dfs2,dfs3]))

In [25]:
## Loading the data:
path = "/home/fatima/Desktop/project_3/"
x = file[:,0]
y = file[:,1]

In [26]:
## plotting the data
plt.style.use('ggplot')
plt.figure(figsize=(8,5))
plt.plot(x,y)
plt.xlabel('Arcseconds')
plt.ylabel('Normalized intensity')
plt.show()

In [27]:
## For the purpose of fitting, we restrict the latter to data points located above x = 0
ind = np.where(x>=0)
x = x[ind]
y = y[ind]
weights = np.sqrt(np.abs(y)) ## Poisson weighting

In [49]:
p0=[0.3, 0.3, 0.2, 1, 2, 3] ## initial guess best-fit parameters
popt, pcov = curve_fit(SL_fit,x,y,p0,method='lm',sigma=weights,jac=DF,ftol=1e-8,xtol=1e-8,maxfev=5000)
chi_sq = np.sum((1/weights**2)*(SL_fit(x,*popt)-y)**2)

In [54]:
print popt #to view the best fit parameters
print chi_sq

[ 0.52750103  0.28882569  0.10191755  0.25905336  0.76540583  2.83343005]
0.00122093770376

In [48]:
xnew = np.linspace(x.min(),x.max(),1000)
y_fit = SL_fit(xnew,*popt)

In [35]:
plt.style.use('ggplot')
plt.figure(figsize=(8,5))
plt.plot(xnew, y_fit,color='blue', label='best fit')
plt.plot(x,y,'r--', label='data')
plt.xlabel('Arcseconds')
plt.ylabel('Normalized intensity')
plt.legend(loc='upper right')
plt.show()


The second method we will use is via scipy.optimize.least_squares . With this module, both he objective function and jacobian matrix have to be defined differently than in scipy.optimize.curve_fit :¶

In [39]:
from scipy.optimize import least_squares

In [40]:
## objective function
def Leastsquares(params, x,y,w):
w1 = params[0]
w2 = params[1]
w3 = params[2]
s1 = params[3]
s2 = params[4]
s3 = params[5]
model = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+(1-w1-w2-w3))
return (1/w)*(model-y)

## as you can see, in the objective function has to return the difference between the modelled data and data points

In [41]:
## jacobian matrix
def Jac(popt,x,y,w):
w1 = popt[0]
w2 = popt[1]
w3 = popt[2]
s1 = popt[3]
s2 = popt[4]
s3 = popt[5]
dfw1 = 0.5*(Erfc(x,s1)-1)
dfw2 = 0.5*(Erfc(x,s2)-1)
dfw3 =  0.5*(Erfc(x,s3)-1)
dfs1 = w1*x*np.exp(-x**2/(2*s1**2))/ (np.sqrt(2*np.pi)*s1**2)
dfs2 = w2*x*np.exp(-x**2/(2*s2**2))/ (np.sqrt(2*np.pi)*s2**2)
dfs3 = w3*x*np.exp(-x**2/(2*s3**2))/ (np.sqrt(2*np.pi)*s3**2)
return np.transpose(np.array(1/w)*([dfw1,dfw2,dfw3,dfs1,dfs2,dfs3]))

In [42]:
## loading the data again:
x = file[:,0]
y = file[:,1]

ind = np.where(x>=0)
x = x[ind]
y = y[ind]
weights = np.sqrt(np.abs(y)) ## Poisson weighting

In [52]:
## calling the objective function and fitting data points
res_leastsquares = least_squares(Leastsquares, p0,method='lm',jac=Jac, args=(x,y,weights),ftol=1e-8,xtol=1e-8,max_nfev=5000)
##computing the chi-square value
chi_sq_2 = np.sum((Leastsquares(res_leastsquares.x,x,y,weights))**2)

In [53]:
## printing best-fit parameters
print res_leastsquares.x
print chi_sq_2

[ 0.52750174  0.28882527  0.10191733  0.2590536   0.76540736  2.83343784]
0.00122093770322


The THIRD method we will use is via scipy.optimize.leasq , which is pretty similar to scipy.optimize.least_squares except for the saved results format¶

In [55]:
from scipy.optimize import leastsq

In [56]:
## objective function
def Leastsq(params, x,y,w):
w1 = params[0]
w2 = params[1]
w3 = params[2]
s1 = params[3]
s2 = params[4]
s3 = params[5]
model = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+(1-w1-w2-w3))
return (1/w)*(model-y)

##jacobian matrix same as in leastsquares method

In [57]:
res_leastsq = leastsq(Leastsq,p0,args=(x,y,weights),Dfun=Jac,ftol=1e-8,xtol=1e-8,maxfev=5000)
chi_sq_3 = np.sum((Leastsq(res_leastsq[0],x,y,weights))**2)

In [58]:
print res_leastsq[0]
print chi_sq_3

[ 0.52750103  0.28882569  0.10191755  0.25905336  0.76540583  2.83343005]
0.00122093770376


The FOURTH method we will use is the powerful package lmfit . LMFIT has the advantage of fixing certain parameters while letting others unchanged, and the possibiliy to introduce bounds on the variables.¶

In [71]:
import lmfit
from lmfit import Minimizer, minimize, Parameters, report_fit

In [63]:
## defining the objective function
def residual(params, x,y,w):
w1 = params['omega1']
w2 = params['omega2']
w3 = params['omega3']
s1 = params['sigma1']
s2 = params['sigma2']
s3 = params['sigma3']

model = 0.5*(w1*Erfc(x,s1)+w2*Erfc(x,s2)+w3*Erfc(x,s3)+(1-w1-w2-w3))

return (1/w)*(y-model)

In [65]:
# using the class Parameters in which we can name the variables and assign first guess values and boundaries, if needed
params = Parameters()

In [67]:
## fitting with lmfit via lmfit.minimize
out = lmfit.minimize(residual,params,method='leastsq',args=(x,y,weights),ftol=1e-8,xtol=1e-8,maxfev=5000)

In [74]:
## printing output
report_fit(out.params)

[[Variables]]
omega1:   0.52750083 +/- 0.003494 (0.66%) (init= 0.3)
omega2:   0.28882580 +/- 0.002845 (0.99%) (init= 0.3)
omega3:   0.10191760 +/- 0.001193 (1.17%) (init= 0.3)
sigma1:   0.25905329 +/- 0.001357 (0.52%) (init= 1)
sigma2:   0.76540540 +/- 0.007152 (0.93%) (init= 2)
sigma3:   2.83342789 +/- 0.047380 (1.67%) (init= 3)
[[Correlations]] (unreported correlations are <  0.100)
C(omega1, sigma2)            =  0.899
C(omega1, omega2)            = -0.895
C(omega1, sigma1)            =  0.889
C(omega3, sigma2)            = -0.870
C(omega3, sigma3)            = -0.853
C(omega2, sigma1)            = -0.851
C(sigma2, sigma3)            =  0.729
C(sigma1, sigma2)            =  0.716
C(omega2, sigma2)            = -0.632
C(omega1, omega3)            = -0.628
C(omega1, sigma3)            =  0.487
C(omega3, sigma1)            = -0.457
C(sigma1, sigma3)            =  0.348
C(omega2, omega3)            =  0.230