[SciPy-user] global fitting of data sets with common parameters
gael.varoquaux at normalesup.org
Mon Sep 25 17:16:17 CDT 2006
> 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*x1[k]+p[1*(2*j+2)])
# Target functions. note the same p for all sets. it is a matrix:
# first index denoting the set, second index denoting a point in a set.
### 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
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
# Legend the plot
title("Oscillations in the compressed trap")
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/p2), fontsize = 16,
transform = ax.transAxes)
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.
More information about the SciPy-user