[SciPy-user] global fitting of data sets with common parameters
Gael Varoquaux
gael.varoquaux at normalesup.org
Mon Sep 25 18:09:03 CDT 2006
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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
More information about the SciPy-user
mailing list