Curve fitting - Bounds

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

## importing the libraries
import numpy as np
import math as m
import scipy
import matplotlib.pyplot as plt
from scipy import special

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

## Importing data
path = "/home/fatima/Desktop/project_3/"
file = np.loadtxt(path+'limb_profile_av_norm_shifted')
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.

from scipy.optimize import curve_fit

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.

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)

print popt

[ 0.52750384  0.28882404  0.10191667  0.25905432  0.76541192  2.83346098]

## 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)

print popt_un

[ 0.52750103  0.28882568  0.10191755  0.25905336  0.76540583  2.83343007]

As you can see, there is no major difference in the output best-fit parameters between the constrained and unconstrained fit. This implies that the solution converges alright to the global one.

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.

from scipy.optimize import least_squares

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)

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

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.

import lmfit
from lmfit import Minimizer, minimize, Parameters, report_fit

## 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)

params = Parameters()
params.add('omega1', value=0.3,min=0,max=1)
params.add('omega2', value=0.3,min=0,max=1)
params.add('omega3', value=0.3,min=0,max=1)
params.add('sigma1', value=0.30,min=0,max=5)
params.add('sigma2', value=2,min=0,max=5)
params.add('sigma3',value=3,min=0,max=5)
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

## for example if we want to fix the parameter omega3 in the params class, we can type the following:
params = Parameters()
params.add('omega1', value=0.3,min=0,max=1)
params.add('omega2', value=0.3,min=0,max=1)
params.add('omega3', value=0.3,min=0,max=1)
params.add('sigma1', value=0.30,vary=False)#,min=0,max=5)
params.add('sigma2', value=2,min=0,max=5)
params.add('sigma3',value=3,min=0,max=5)
out2 = lmfit.minimize(residual,params,method='leastsq',args=(x,y,weights))
out2.params.pretty_print()

Name       Value      Min      Max   Stderr     Vary     Expr
omega1    0.6164        0        1 0.002565     True     None
omega2    0.2268        0        1 0.003173     True     None
omega3   0.08511        0        1 0.001487     True     None
sigma1       0.3     -inf      inf        0    False     None
sigma2     0.973        0        5  0.01644     True     None
sigma3     3.985        0        5   0.3694     True     None