analysis: gauss_deriv_fit.py Source File

00001 "Fit the derivative of a Gaussian to a data set.  Mostly an example of using fitting_toolkit.py"
00002 _rcsid="$Id: gauss_deriv_fit.py,v 1.3 2003/05/30 13:31:55 mendenhall Release-20030716 $"
00003 
00004 import fitting_toolkit
00005 import Numeric
00006 from math import sqrt
00007 
00008 class gauss_deriv_fit(fitting_toolkit.fit):
00009     "fit a constant baseline + (x-mu)*gaussian y=y0+a*(x-mu)*exp( -(x-xmu)**2/(2*xsig**2) )"
00010     def function(self, p, r):
00011         z0, a, xmu, xsigma = p
00012         xsigi=-1.0/(2.0*xsigma**2)
00013         return z0+a*(r-xmu)*Numeric.exp(xsigi*(r-xmu)**2)
00014 
00015     def derivs(self): 
00016         #analytic derivatives for a 1-d gaussian*(x-mu)
00017         #z0+a*(x-mu)*exp( -(x-xmu)**2/(2*xsig**2) )
00018         z0, a, xmu, xsigma = self.funcparams
00019         n=self.pointcount
00020         x=self.xarray[0,:n]
00021         xsigi=-1.0/(2.0*xsigma**2)
00022         dx=x-xmu
00023         dx2=dx*dx
00024         expfact=Numeric.exp(xsigi*dx2)
00025         z=a*expfact*dx
00026         
00027         dd = Numeric.zeros((n, 4), self.atype)
00028         dd[:,0]=1.0
00029         dd[:,1]=expfact*dx
00030         dd[:,2]=(-2.0*xsigi*dx*dx - 1)*a*expfact
00031         dd[:,3]=(-2.0*xsigi/xsigma)*(dx2*z)
00032         
00033         return dd   
00034 
00035 if __name__=="__main__":
00036     
00037     x=gauss_deriv_fit()
00038 
00039     z0, a, xmu, xsigma =  1., 15., 73., 10.
00040     x.set_initial_params([z0+5, a-10, xmu+5, xsigma-5])
00041     
00042     xlist=Numeric.array(range(100),Numeric.Float)
00043     ylist=x.compute_funcvals(params=[z0, a, xmu, xsigma], xvals=xlist)
00044     
00045     x.add_points(xlist, ylist)
00046     
00047     print "\n\n***Start nonlinear test fit***" 
00048     for i in range(10):
00049         x.lm_fit_step()
00050         print Numeric.array_str(x.funcparams, precision=5),  sqrt(x.reduced_chi2)
00051