[SciPy-user] using errors in scipy.optimize.leastsq

Tommy Grav tgrav@mac....
Sat Apr 14 17:47:56 CDT 2007


I am using scipy.optimize.leastsq to fit a sinusoid function to a set of
observations given by three arrays atime,amag,aerr. I am able to make
the leastsq function find what seems to be the appropriate solution,
but I am now wondering how to use the amerr (the array of errors in
the measurement to weight the individual points used in the solution.
Does anyone have any hints for doing this?

Cheers
   Tommy

Snipped example:

def func(x,period,time,mag,merr,n):
     diff = ndarray(len(time))
     tmp = 2.*pi/period
     for l in range(1,n+1):
         fit = x[0] + x[2]*sin(tmp*(time - x[1])) + x[3]*cos(tmp* 
(time - x[1]))
     diff = fit - mag
     return diff


mavg = average(amag)
mmax = max(amag) - mavg
mmin = mavg - min(amag)
print mavg, mmax, mmin

x0 = array([mavg,54128.0,mmax,mmin])
norder = 1
plist = array([float(p/1000.) for p in range(10,1000)])
chilist = ndarray(len(plist))
n = 0
for period in plist:
     res = leastsq(func,x0,args= 
(period,atime,amag,amerr,norder),full_output=1,col_deriv=1,epsfcn=0.001)
     diff = func(res[0],period,atime,amag,amerr,norder)
     chilist[n] = sqrt(dot(diff,diff)/(len(diff)-1.))
     n += 1


period = plist[chilist.argmin()]
print "Period = %12.5f h" % (period*24.)

res = leastsq(func,x0,args= 
(period,atime,amag,amerr,norder),full_output=1,col_deriv=1)
print "Amplitude = %12.5f mag" % (res[0][2])



More information about the SciPy-user mailing list