global fitting of data sets with common parameters
Hello I have data sets and would like to fit them globally with common parameters. Is there a way to do this using existing scipy modules leastsq. Any pointers for algo and suggestions are welcome. Thanks Morovia.
On Wed, Aug 09, 2006 at 10:15:04AM +0200, morovia morovia wrote:
I have data sets and would like to fit them globally with common parameters. Is there a way to do this using existing scipy modules leastsq.
Any pointers for algo and suggestions are welcome.
Hi, I just uploaded a example to do this in the cookbook, on scipy.org . Does it suit your needs ? Gaël
Gael Varoquaux wrote:
On Wed, Aug 09, 2006 at 10:15:04AM +0200, morovia morovia wrote:
... fit ... globally with common parameters. ... I just uploaded a example to do this in the cookbook, on scipy.org . ...
Hey Gaël! That's absolutely great! I was about to write such a thing myself. I guess you simply add up the chi^2 of all data sets and get this minimized? Or does one have to use the reduced chi^2? Could you just please specify a bit more in detail /where/ you uploaded your example in the cookbook? Sorry for being that blind, Sebastian.
On Mon, Sep 25, 2006 at 10:39:51AM +0200, Sebastian Busch wrote:
That's absolutely great! I was about to write such a thing myself. I guess you simply add up the chi^2 of all data sets and get this minimized? Or does one have to use the reduced chi^2?
Could you just please specify a bit more in detail /where/ you uploaded your example in the cookbook?
http://scipy.org/Cookbook/FittingData , linked on the main cookbook page, near the top (I thought this was a really useful example). If you think you can add chi^2 without rendering the figure unreadable, please do. Gaël
Gael Varoquaux wrote:
On Mon, Sep 25, 2006 at 10:39:51AM +0200, Sebastian Busch wrote:
... please specify ... /where/ you uploaded your example ...
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()
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
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
Gael Varoquaux wrote:
Gael Varoquaux wrote:
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.
hey there! Thanks a lot! This is great! The only thing you might want to change (I am not allowed to) is the leading space in line 18. I am not sure what does not seem good to you -- the for-loops etc (I totally agree) or the common fitting (in this case, I would be happy if you would like to discuss this in more detail). Thanks again, Sebastian.
participants (3)
-
Gael Varoquaux -
morovia morovia -
Sebastian Busch