[SciPy-user] global fitting of data sets with common parameters

Sebastian Busch sebastian_busch at gmx.net
Mon Sep 25 15:16:19 CDT 2006


Gael Varoquaux wrote:
> On Mon, Sep 25, 2006 at 10:39:51AM +0200, Sebastian Busch wrote:
>> ... please specify ... /where/ you uploaded your example ...
> 
> http://scipy.org/Cookbook/FittingData ...

Hi Gaël!

Thanks for your reply. I think this is not quite what the OP wanted. In
any case, it is not what I am looking for ;)

You fit two data sets with two different functions. Completely
independent. You could do the same in two runs of your program working
on the different data sets.

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.

So you sit down and want to fit the data. Not in two runs, then
averaging the two frequencies you get, but in one big fit. This is
better because one measurement was more precise than the other and it's
hard to account for that after the fit finished.

I tried to modify your example in that way. As I am pretty much a bloody
beginner in python, I failed. I post my try anyway to show you more
clearly what I am thinking of. And, I admit, in the hope you can do
better ;)

Thanks for your time,
Sebastian.

+++++++++++++++++++

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...


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





# 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

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()



More information about the SciPy-user mailing list