[SciPy-User] curve_fit and least squares

Warren Weckesser warren.weckesser@enthought....
Thu Oct 22 15:12:10 CDT 2009


Kris,

Your script worked for me if I explicitly converted everything to numpy 
arrays.  Here's my edited version:

----------
#!/usr/bin/env python
import numpy as np
import scipy as sp
import pylab as pl
from scipy.optimize.minpack import curve_fit

x = np.array([  50.,  110.,  170.,  230.,  290.,  350.,  410.,  470.,  
530.,  590.])
y = np.array([ 3173.,  2391.,  1726.,  1388.,  1057.,   786.,   598.,   
443.,   339.,   263.])

smoothx = np.linspace(x[0], x[-1], 20)
guess_a, guess_b, guess_c = 4000, -0.005, 100
guess = [guess_a, guess_b, guess_c]

f_theory1 = lambda t, a, b, c: a * np.exp(b * t) + c

p, cov = curve_fit(f_theory1, x, y, p0=np.array(guess))

pl.clf()
f_fit1 = lambda t: p[0] * np.exp(p[1] * t) + p[2]
# pl.plot(x, y, 'b.', smoothx, f_theory1(smoothx, guess_a, guess_b, 
guess_c))
pl.plot(x, y, 'b.', smoothx, f_fit1(smoothx), 'r-')
pl.show()

##
## EOF
##

----------

Warren


Kris Maynard wrote:
> Hi,
>
> Thanks for your responses. After some more digging and some more 
> testing I'm beginning to think that the algorithm used by curve_fit 
> simply isn't robust enough for the data that I am trying to fit. Below 
> is an example of some experimental radioactive decay data that I am 
> trying to fit to an exponential decay. 
>
> #!/usr/bin/env python
> import numpy as np
> import scipy as sp
> import pylab as pl
> from scipy.optimize.minpack import curve_fit
>
> x = [  50.,  110.,  170.,  230.,  290.,  350.,  410.,  470.,  530.,  590.]
> y = [ 3173.,  2391.,  1726.,  1388.,  1057.,   786.,   598.,   443.,   
> 339.,   263.]
>
> smoothx = np.linspace(x[0], x[-1], 20)
> guess_a, guess_b, guess_c = 4000, -0.005, 100
> guess = [guess_a, guess_b, guess_c]
>
> f_theory1 = lambda t, a, b, c: a * np.exp(b * t) + c
> f_theory2 = lambda t, a, b: np.exp(a * t) + b
>
> pl.plot(x, y, 'b.', smoothx, f_theory1(smoothx, guess_a, guess_b, 
> guess_c))
> pl.show()
>
> p, cov = curve_fit(f_theory1, x, y)
> #p, cov = curve_fit(f_theory2, x, y)
>
> # the following gives:
> #   ValueError: shape mismatch: objects cannot be broadcast to a 
> single shape
> #p, cov = curve_fit(f_theory1, x, y, p0=guess)
>
> pl.clf()
> f_fit1 = lambda t: p[0] * np.exp(p[1] * t) + p[2]
> #f_fit2 = lambda t: np.exp(p[0] * t) + p[1]
> pl.plot(x, y, 'b.', smoothx, f_fit1(smoothx), 'k-')
> pl.show()
>
> ##
> ## EOF
> ##
>
> As you can see, I have tried to fit using 2 or 3 parameters with no 
> luck. Is there something that I could do to make this work? I have 
> tried this exact thing in matlab and it worked the first time. 
> Unfortunately, I would really like to use python as I find it in 
> general more intuitive than matlab.
>
> Thanks,
> ~Kris
>
> On Wed, Oct 7, 2009 at 9:40 AM, Bruce Southey <bsouthey@gmail.com 
> <mailto:bsouthey@gmail.com>> wrote:
>
>     On Wed, Oct 7, 2009 at 1:19 AM,  <josef.pktd@gmail.com
>     <mailto:josef.pktd@gmail.com>> wrote:
>     > On Wed, Oct 7, 2009 at 1:36 AM, Kris Maynard <maynard@bu.edu
>     <mailto:maynard@bu.edu>> wrote:
>     >> I am having trouble with fitting data to an exponential curve.
>     I have an x-y
>     >> data series that I would like to fit to an exponential using
>     least squares
>     >> and have access to the covariance matrix of the result. I
>     summarize my
>     >> problem in the following example:
>     >>
>     >> import numpy as np
>     >> import scipy as sp
>     >> from scipy.optimize.minpack import curve_fit
>     >>
>     >> A, B = 5, 0.5
>     >> x = np.linspace(0, 5, 10)
>     >> real_f = lambda x: A * np.exp(-1.0 * B * x)
>     >> y = real_f(x)
>     >> ynoisy = y + 0.01 * np.random.randn(len(x))
>     >>
>     >> exp_f = lambda x, a, b: a * np.exp(-1.0 * b * x)
>     >>
>     >> # this line raises the error:
>     >>
>     >> #  RuntimeError: Optimal parameters not found: Both
>     >>
>     >> #  actual and predicted relative reductions in the sum of squares
>     >>
>     >> #  are at most 0.000000 and the relative error between two
>     >>
>     >> #  consecutive iterates is at most 0.000000
>     >>
>     >> params, cov = curve_fit(exp_f, x, ynoisy)
>     >
>
>     Could you please first plot your data?
>     As you would see, the curve is very poorly defined with those model
>     parameters and range. So you are asking a lot from your model and
>     data. At least you need a wider range with those parameters or Josef
>     says different parameter(s):
>
>     > this might be the same as
>      http://projects.scipy.org/scipy/ticket/984 and
>     > http://mail.scipy.org/pipermail/scipy-user/2009-August/022090.html
>     >
>     > If I increase your noise standard deviation from 0.1 to 0.2 then
>     I do get
>     > correct estimation results in your example.
>     >
>     >>
>     >> I have tried to use the minpack.leastsq function directly with
>     similar
>     >> results. I also tried taking the log and fitting to a line with
>     no success.
>     >> The results are the same using scipy 0.7.1 as well as
>     0.8.0.dev5953. Am I
>     >> not using the curve_fit function correctly?
>     >
>     > With   minpack.leastsq   error code 2 should be just a warning.
>     If you get
>     > incorrect parameter estimates with optimize.leastsq, besides the
>     warning, could
>     > you post the example so I can have a look.
>     >
>     > It looks like if you take logs then you would have a problem
>     that is linear in
>     > (transformed) parameters, where you could use linear least
>     squares if you
>     > just want a fit without the standard errors of the original
>     parameters
>     > (constant)
>
>     The errors will be multiplicative rather than additive.
>
>     Bruce
>
>     >
>     > I hope that helps.
>     >
>     > Josef
>     >
>     >
>     >> Thanks,
>     >> ~Kris
>     >> --
>     >> Heisenberg went for a drive and got stopped by a traffic cop.
>     The cop asked,
>     >> "Do you know how fast you were going?" Heisenberg replied, "No,
>     but I know
>     >> where I am."
>     >>
>     >> _______________________________________________
>     >> SciPy-User mailing list
>     >> SciPy-User@scipy.org <mailto:SciPy-User@scipy.org>
>     >> http://mail.scipy.org/mailman/listinfo/scipy-user
>     >>
>     >>
>     > _______________________________________________
>     > SciPy-User mailing list
>     > SciPy-User@scipy.org <mailto:SciPy-User@scipy.org>
>     > http://mail.scipy.org/mailman/listinfo/scipy-user
>     >
>
>     Hi,
>     _______________________________________________
>     SciPy-User mailing list
>     SciPy-User@scipy.org <mailto:SciPy-User@scipy.org>
>     http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
>
> -- 
> Heisenberg went for a drive and got stopped by a traffic cop. The cop 
> asked, "Do you know how fast you were going?" Heisenberg replied, "No, 
> but I know where I am."
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>   



More information about the SciPy-User mailing list