OK, I completed to cookbook page. This may not be the most efficient way to do things, but I find it quite readable. GaÃl On Tue, Sep 26, 2006 at 12:16:17AM +0200, Gael Varoquaux wrote:
What makes your example similar to what I want is the frequency: it is nearly the same. Now let's assume these are two motions of which you know that the frequencies should be the same. Any differences are due to measurement errors.
OK, I'll first comment your code (I shaped it to be able to read it better. Your line are to long, especially in comments, and your comments should be more often on separate lines, indented the same way as the code). My comments are prefixed by ###
+++++++++++++++++ from pylab import * from scipy import *
# Generate data points with noise num_points = 150
# as i got confused with xT, Xt, ... it's now data set 1 and 2, each set # having x/y coordinates. x1 = linspace(5., 8., num_points) x2 = x1 xexp = array([x1, x2]) # i put them together. not sure whether this is the right way... ### You probably want to you use c_[x1,x2]
y1 = 11.86*cos(2*pi/0.81*x1-1.32) + 6.*rand(num_points) y2 = -32.14*cos(2*pi/0.8*x2-1.94) + 3.*rand(num_points) yexp = [y1, y2] # same thing
### !!!! I thing you misstyped here you meant yexp = c_[y1,y2]
# Fit setup def ytheo (p,x) : # yth = zeros((x1,x1)) how to initialize? for j in range(len(xexp)): # count through data sets for k in range(len(x1)): # count through the data points within a (the first) set yth[j,k] = p[(2*j+1)]*cos(2*pi/p[0]*x1[k]+p[1*(2*j+2)]) # Target functions. note the same p[0] for all sets. it is a matrix: # first index denoting the set, second index denoting a point in a set. return yth
### A for loop !! Your are kidding. You should try to avoid this as much ### as possible, and work only on arrays. For loops cost a lot of time.
errfunc = lambda p, x, y: ytheo(p,x) - yexp # Distance to the target function
# fit p0 = c_[0.8, -15., 0., -15., 0.] # Initial guess for the parameters # 0: common "frequency"; 1: 1st set ampl; 2: 1st set phase; # 3: 2nd set ampl; 4: 2nd set phase
p1,success = optimize.leastsq(errfunc, p0.copy(), args = (xexp, yexp))
time = linspace(x1.min(),x1.max(),100) plot(x1,y1,"ro",time,fitfunc(p1,time),"r-") # Plot of the data and the fit plot(x2,y2,"b^",time,fitfunc(p2,time),"b-")
# Legend the plot title("Oscillations in the compressed trap") xlabel("time [ms]") ylabel("displacement [um]") legend(('1 position', '1 fit', '2 position', '2 fit'))
ax = axes()
text(0.8, 0.07,'x freq : %.3f kHz \n y freq : %.3f kHz' %(1/p1[1],1/p2[1]), fontsize = 16, horizontalalignment='center', verticalalignment='center', transform = ax.transAxes)
show() +++++++++++++++++
I'll upload an example show you how I would what you are trying to do. The wiki is just so much more readable than mail :->. I don't exactly see what your are trying to do, but it does not seems to good to me.
GaÃl
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user