[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