[SciPy-User] curve_fit and least squares
Kris Maynard
maynard@bu....
Thu Oct 22 14:54:14 CDT 2009
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> wrote:
> On Wed, Oct 7, 2009 at 1:19 AM, <josef.pktd@gmail.com> wrote:
> > On Wed, Oct 7, 2009 at 1:36 AM, Kris Maynard <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/984and
> > 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
> >> http://mail.scipy.org/mailman/listinfo/scipy-user
> >>
> >>
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User@scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
>
> Hi,
> _______________________________________________
> SciPy-User mailing list
> 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."
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20091022/4167d281/attachment.html
More information about the SciPy-User
mailing list