Python Tutorial 3

import numpy as np
import scipy as sp
import pylab as plt
import scipy.stats as stats
get_ipython().magic(u'matplotlib inline')


# ###Nonlinear fitting
#
# Here we'll use the scipy.optimize.curve_fit function, which implements the L-M nonlinear fitting method, to fit a  model (a gaussian with linear polynomial) to some fake data.

# In[36]:

def model(x,a,b,c,d,e):
    return a+b*x+c*np.exp(-(x-d)**2/2/e**2)

#Make fake data
x = np.arange(-5,5,.5)
err = np.ones(len(x))*1.2
y = np.random.normal(model(x,10,2,20,1,1),scale=err)
plt.errorbar(x,y,err,fmt='o')

#Fit the data
p,cov = sp.optimize.curve_fit(model,x,y,sigma=err,absolute_sigma=True,p0=[5,1,25,1.5,2])
for i in range(len(p)):
    print p[i],'+/-',np.sqrt(cov[i][i])
plt.plot(x,model(x,*p))
plt.plot(x,p[0]+p[1]*x)
plt.plot(x,p[2]*np.exp(-(x-p[3])**2/2/p[4]**2))
plt.show()

chisq = np.sum((y-model(x,*p))**2/err**2)
dof = len(y)-len(p)
print 'chisq =', chisq
print 'dof', dof
print 'PTE', stats.chisqprob(chisq,dof)