[SciPy-User] curve_fit and least squares

Skipper Seabold jsseabold at gmail.com
Thu Oct 22 16:09:04 EDT 2009 

On Thu, Oct 22, 2009 at 3:54 PM, Kris Maynard <maynard at bu.edu> 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


I didn't look too closely, but the error message suggests that your
fit_theory functions return an array that is not the right dimension
compared to y

Consider:

f_theory1(smoothx, guess_a, guess_b, guess_c) - y

#ValueError: shape mismatch: objects cannot be broadcast to a single shape

f_theory1(smoothx, guess_a, guess_b, guess_c).shape
# (20,)
# of type array

type(y)
# list

This return an array, but your y is still a list, so

In [52]: curve_fit(f_theory1, np.array(x), np.array(y), p0=(guess_a,
guess_b, guess_c))
Out[52]:
(array([  3.97627120e+03,  -4.76859072e-03,   2.93930899e+01]),
 array([[  2.03195757e+03,  -1.17944031e-03,  -3.33962266e+02],
       [ -1.17944031e-03,   3.14507326e-08,  -6.91762164e-03],
       [ -3.33962266e+02,  -6.91762164e-03,   1.80377505e+03]]))

I have no idea if this result makes sense, but it might point you in
the right direction.

Skipper





> On Wed, Oct 7, 2009 at 9:40 AM, Bruce Southey <bsouthey at gmail.com> wrote:
>>
>> On Wed, Oct 7, 2009 at 1:19 AM,  <josef.pktd at gmail.com> wrote:
>> > On Wed, Oct 7, 2009 at 1:36 AM, Kris Maynard <maynard at 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 at scipy.org
>> >> http://mail.scipy.org/mailman/listinfo/scipy-user
>> >>
>> >>
>> > _______________________________________________
>> > SciPy-User mailing list
>> > SciPy-User at scipy.org
>> > http://mail.scipy.org/mailman/listinfo/scipy-user
>> >
>>
>> Hi,
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>

More information about the SciPy-User mailing list