Hi, Sorry for the bunch of questions, but I'm preparing some 'hands-on' on scipy, and I'm excercising many parts of it. I've tried to follow the Travis' tutorial for finding the best fit of a series of data, but my results are different from those stated there: In [96]: x = arange(0,6e-2,6e-2/30) In [97]: A,k,theta = 10, 1.0/3e-2, pi/6 In [98]: y_true = A*sin(2*pi*k*x+theta) In [99]: y_meas = y_true + 2*randn(len(x)) In [100]: def residuals(p, y, x): .....: A,k,theta = p .....: err = y-A*sin(2*pi*k*x+theta) .....: return err .....: In [101]: p0 = [8, 1/2.3e-2, pi/3] In [107]: optimize.leastsq(residuals, p0, args=(y_meas, x)) Out[107]: (array([ -2.17634575, 57.4116454 , -3.03266312]), 'Both actual and predicted relative reductions in the sum of squares\n are at most 0.000000') this results are quite different from the 'true' values: [ 10. 33.3333 0.5236] Do you have any idea about this error? I suppose that the message should indicate it, but I can't understand it. I'm using scipy 0.3.2 on Linux. BTW, thanks for clearifying the question about derivative() parameters. It's much more clear now. -- Francesc Alted
Francesc Alted wrote:
Hi,
Sorry for the bunch of questions, but I'm preparing some 'hands-on' on scipy, and I'm excercising many parts of it.
I've tried to follow the Travis' tutorial for finding the best fit of a series of data, but my results are different from those stated there:
In [96]: x = arange(0,6e-2,6e-2/30)
In [97]: A,k,theta = 10, 1.0/3e-2, pi/6
In [98]: y_true = A*sin(2*pi*k*x+theta)
In [99]: y_meas = y_true + 2*randn(len(x))
In [100]: def residuals(p, y, x): .....: A,k,theta = p .....: err = y-A*sin(2*pi*k*x+theta) .....: return err .....:
In [101]: p0 = [8, 1/2.3e-2, pi/3]
In [107]: optimize.leastsq(residuals, p0, args=(y_meas, x)) Out[107]: (array([ -2.17634575, 57.4116454 , -3.03266312]), 'Both actual and predicted relative reductions in the sum of squares\n are at most 0.000000')
this results are quite different from the 'true' values: [ 10. 33.3333 0.5236]
Do you have any idea about this error? I suppose that the message should indicate it, but I can't understand it. I'm using scipy 0.3.2 on Linux.
Can you look at the output of peval(x,p) and compare it to y_true? I just re-ran this experiment and get array([-10.053690336957587, 33.659758123232756, 3.591816659730455)]) which seems right (except for the fact that some of the phase has been used to give -10 instead of 10). -Travis O.
A Dijous 21 Octubre 2004 21:01, Travis Oliphant va escriure:
Can you look at the output of peval(x,p) and compare it to y_true?
I just re-ran this experiment and get array([-10.053690336957587, 33.659758123232756, 3.591816659730455)])
which seems right (except for the fact that some of the phase has been used to give -10 instead of 10).
Mmmm, I've re-run the example again and get [-10.22083869, 34.03240033, 3.41523692] so, I'm afraid that something weird happened in the run I've made yesterday. Perhaps (I'm guessing here) the way to compute deviations, i.e.: y_meas = y_true + 2*randn(len(x)) can give place to values too far from y_true and that can make the algorithm to fail. It's a pity, but I dind't saved the values from my experiment yesterday, so I can't reproduce the problem. Anyway, maybe it would be sensible to change the line where the computations of deviations are made to read: y_meas = y_true + randn(len(x)) although this may not be the cause, and it was just a human error :-/ Thanks, -- Francesc Alted
Hi,
Sorry for the bunch of questions, but I'm preparing some 'hands-on' on scipy, and I'm excercising many parts of it.
Here is a different simple scipy least squares fit: http://www.scipy.org/mailinglists/mailman?fn=scipy-user/2004-January/002452.... James Phillips http://zunzun.com
Hi James, A Dijous 21 Octubre 2004 22:12, James R. Phillips va escriure:
Here is a different simple scipy least squares fit:
http://www.scipy.org/mailinglists/mailman?fn=scipy-user/2004-January/002452....
Thanks for your suggestion. Here is another version of your code, just a bit more 'vectorized': ----------------------------------------------------------------------- from scipy import * xdata = array([5.357, 5.457, 5.797, 5.936, 6.161, 6.697, 6.731, 6.775, 8.442, 9.769, 9.861]) ydata = array([0.376, 0.489, 0.874, 1.049, 1.327, 2.054, 2.077, 2.138, 4.744, 7.068, 7.104]) transpose(array([[1]*11, xdata, xdata*xdata])) # y = 1 + x + x*x coeffs = linalg.basic.lstsq(matrix, ydata)[0] print "scipy.linalg.basic.lstsq curve fitting example" print "fitting data to quadratic equation y = a + bx + cx^2" ycalc = coeffs[0] + coeffs[1] * xdata + coeffs[2] * xdata * xdata error = ycalc - ydata print "yields: x data y data calc value error" for i in range(len(xdata)): print " % .3f % .3f % .3f % .3f" % (xdata[i], ydata[i], ycalc[i], error[i]) print ------------------------------------------------------------------------- Regards, -- Francesc Alted
Thanks for your suggestion. Here is another version of your code, just a bit more 'vectorized':
transpose(array([[1]*11, xdata, xdata*xdata])) # y = 1 + x + x*x
The comment at the end of this line is not quite correct, I think. Look at the one in the original example again and see what you make of the difference. James Phillips http://zunzun.com
Ooops, a small glitch went in my 'vectorized' version. This has been tested: from scipy import * xdata = array([5.357, 5.457, 5.797, 5.936, 6.161, 6.697, 6.731, 6.775, 8.442, 9.769, 9.861]) ydata = array([0.376, 0.489, 0.874, 1.049, 1.327, 2.054, 2.077, 2.138, 4.744, 7.068, 7.104]) matrix=transpose(array([[1]*11, xdata, xdata*xdata])) # y = 1 + x + x*x coeffs = linalg.basic.lstsq(matrix, ydata)[0] print "scipy.linalg.basic.lstsq curve fitting example" print "fitting data to quadratic equation y = a + bx + cx^2" ycalc = coeffs[0] + coeffs[1] * xdata + coeffs[2] * xdata * xdata error = ycalc - ydata print "yields: x data y data calc value error" for i in range(len(xdata)): print " % .3f % .3f % .3f % .3f" % (xdata[i], ydata[i], ycalc[i], error[i]) print -- Francesc Alted
participants (3)
-
Francesc Alted -
James R. Phillips -
Travis Oliphant