[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