[SciPy-User] curve_fit and least squares
josef.pktd at gmail.com
josef.pktd at gmail.com
Thu Oct 22 16:29:09 EDT 2009
On Thu, Oct 22, 2009 at 4:20 PM, Kris Maynard <maynard at bu.edu> wrote:
> 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."
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
I haven't looked, but from your discussion, it looks like there is an
y = np.asarray(y) missing in curve_fit ?
Josef
More information about the SciPy-User
mailing list