[SciPy-User] curve_fit and least squares
Kris Maynard
maynard at bu.edu
Thu Oct 22 16:20:50 EDT 2009
Oops.
Right you both are. I suppose the answer to my general confusion is that
jimmying the initial fit parameters in the right way will make curve_fit
work. Thanks!
~Kris
On Thu, Oct 22, 2009 at 4:12 PM, Warren Weckesser <
warren.weckesser at enthought.com> wrote:
> 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 at gmail.com
> > <mailto:bsouthey at gmail.com>> wrote:
> >
> > On Wed, Oct 7, 2009 at 1:19 AM, <josef.pktd at gmail.com
> > <mailto:josef.pktd at gmail.com>> wrote:
> > > On Wed, Oct 7, 2009 at 1:36 AM, Kris Maynard <maynard at bu.edu
> > <mailto: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 <mailto:SciPy-User at scipy.org>
> > >> http://mail.scipy.org/mailman/listinfo/scipy-user
> > >>
> > >>
> > > _______________________________________________
> > > SciPy-User mailing list
> > > SciPy-User at scipy.org <mailto:SciPy-User at scipy.org>
> > > http://mail.scipy.org/mailman/listinfo/scipy-user
> > >
> >
> > Hi,
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User at scipy.org <mailto: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
> >
>
> _______________________________________________
> 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."
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20091022/cb3321bb/attachment.html>
More information about the SciPy-User
mailing list