{[ promptMessage ]}

Bookmark it

{[ promptMessage ]}

lec13_nonlinear_fitting

# lec13_nonlinear_fitting - chisqr...

This preview shows pages 1–2. Sign up to view the full content.

from pylab import * from scipy.optimize import curve_fit def avg(x): return sum(x)/float(len(x)) def sst(x): sum=0.0 xbar=avg(x) for i in range(len(x)): sum += (x[i]-xbar)**2 return sum def ssr(x,y,model,params): sum = 0.0 for i in range(len(y)): sum += (y[i] - model(x[i],*params))**2 return sum def funct(x,a,k,tau,phi): return a*exp(-x/tau)*sin(k*x+phi) # errs=0.4 #Generated simulated data par=[2.,3.0,4.0,0.5] xdata=arange(0,10,0.2) ydata=funct(xdata,*par)+errs*randn(len(xdata)) #initial guess for parameters par0=(10.5,100.0,1.0,0.1) #Perform fit fit=curve_fit(funct,xdata,ydata,p0=par0) #Parse output and extract some statistics fitparams=fit[0] a_fit=fitparams[0] k_fit=fitparams[1] tau_fit=fitparams[2] phi_fit=fitparams[3] #Fit statistics cov=fit[1] #Covariance matrix dof=len(xdata)-len(fitparams) #degrees of freedom chisqr=ssr(xdata,ydata,funct,fitparams) #chi-square is sum of squares of diagonal covariance errbars=20 #Computer uncertainties in estimated parameters from covariance matrix and reduced

This preview has intentionally blurred sections. Sign up to view the full version.

View Full Document
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: chisqr del_A=sqrt(cov[0,0]*sqrt(chisqr/dof)) del_omega=sqrt(cov[1,1]*sqrt(chisqr/dof)) del_tau=sqrt(cov[2,2]*sqrt(chisqr/dof)) del_phi=sqrt(cov[3,3]*sqrt(chisqr/dof)) #Error in parameter p[i] is then: cov[i,i]*sqrt(chisqr/dog) #Generate fit curve xfit=linspace(min(xdata),max(xdata),100) yfit=funct(xfit,a_fit,k_fit,tau_fit,phi_fit) #Plot it all plot(xdata,ydata,'o',label='Data') plot(xfit,yfit,'-',label='Fit') xlabel('Time (s)') ylabel('Amplitude (cm)') legend() print 'Fit Results (error)' print '+'*30 print 'A= %2.3f +/- %1.3f cm' % (a_fit, del_A) print 'omega= %2.3f +/- %1.3f Hz' % (k_fit, del_omega) print 'tau= %2.3f +/- %1.3f sec' % (tau_fit, del_tau) print 'phi= %2.3f +/- %1.3f radians' % (phi_fit, del_phi) r=sqrt(1.0-ssr(xdata,ydata,funct,fitparams)/sst(ydata)) rsq=r**2 print '+'*30 print 'Correlation Coefficient (R^2): %2.3f' % rsq...
View Full Document

• Spring '09