# Model fitting and introducing bounds

## Curve fitting - Bounds¶

#### In this second part, I will explain how to fit a set of data points but with the use of bounds.¶

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

In [3]:
## Importing data
path = "/home/fatima/Desktop/project_3/"
x = file[:,0]
y = file[:,1]
ind = np.where(x>=0)
x = x[ind]
y = y[ind]
weights = np.sqrt(np.abs(y))


#### The first option is to use scipy.optimize.curve_fit. The defined bounds should be in 2 tuples of arrays. The first array should include the lower boundaries of the fit parameters while the second array should include the maximum boundaries.¶

In [4]:
from scipy.optimize import curve_fit

In [5]:
bounds = ((0,0,0,0,0,0),(1,1,1,5,5,5))


#### usually for constrained fitting, the method to be passed to curve_fit is 'trf' instead of 'lm', which is only suitable for unconstrained fit.¶

In [7]:
p0=[0.3, 0.3, 0.2, 1, 2, 3]
popt, pcov = curve_fit(SL_fit,x,y,p0,bounds=bounds,method='trf',sigma=weights,ftol=1e-8,xtol=1e-8,maxfev=5000)

In [8]:
print popt

[ 0.52750384  0.28882404  0.10191667  0.25905432  0.76541192  2.83346098]

In [9]:
## unconstrained fit:
popt_un, pcov_un = curve_fit(SL_fit,x,y,p0,method='lm',sigma=weights,ftol=1e-8,xtol=1e-8,maxfev=5000)

In [10]:
print popt_un

[ 0.52750103  0.28882568  0.10191755  0.25905336  0.76540583  2.83343007]


#### The second option is to use scipy.optimize.least_squares. But then as we saw in curve_fitting we have to redefine the objective function.¶

In [12]:
from scipy.optimize import least_squares

In [15]:
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)

In [17]:
res_leastsquares = least_squares(Leastsquares, p0,method='trf',bounds=bounds,args=(x,y,weights),ftol=1e-8,xtol=1e-8,max_nfev=5000)

In [18]:
print res_leastsquares.x

[ 0.52750384  0.28882404  0.10191667  0.25905432  0.76541192  2.83346098]


#### The third option is to use lmfit . The bounds could be entered as minimum and maximum parameters of the Parameters object.¶

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

In [20]:
## 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 [22]:
params = Parameters()
out = lmfit.minimize(residual,params,method='leastsq',args=(x,y,weights))
#report_fit(out.params)
out.params.pretty_print()

Name       Value      Min      Max   Stderr     Vary     Expr
omega1    0.5275        0        1 0.003494     True     None
omega2    0.2888        0        1 0.002846     True     None
omega3    0.1019        0        1 0.001194     True     None
sigma1    0.2591        0        5 0.001357     True     None
sigma2    0.7654        0        5 0.007153     True     None
sigma3     2.833        0        5  0.04738     True     None


#### One cool application of lmfit is that it allows to fix one or more parameters while varying the others. This can be done by using the parameter vary which set by default to True¶

In [24]:
## for example if we want to fix the parameter omega3 in the params class, we can type the following:
params = Parameters()

Name       Value      Min      Max   Stderr     Vary     Expr