[SciPy-User] fitting cosine with curve_fit - problem with frequency

josef.pktd@gmai... josef.pktd@gmai...
Wed Oct 12 19:55:17 CDT 2011


On Wed, Oct 12, 2011 at 8:23 PM, Brian Blais <bblais@bryant.edu> wrote:
> Hello,
>
> So I am trying to fit some data that is a mixture of oscillations, so I figured I'd try to do a "toy" problem...and that problem is stumping me (full code at the bottom).  I generate data like:
>
> x=pylab.linspace(0,50,100)
> y=cos(x/2)+pylab.randn(len(x))*.1
>
> and fit a model like:
>
> def m1(t,a,b,c):
>    return a*cos(b*t+c)
>
>
> with scipy.optimize.curve_fit.  It seems as if the method used in curve_fit doesn't like oscillatory data, especially when estimating the frequency.  Unless my initial guess for the frequency is *very* close to the right answer (i.e. I need the right answer already to get the answer), it doesn't even get close.
>
> Is there a better way to fit this sort of a function?  Should I do an fft, pick off frequencies, and use those as the initial estimates?  Am I doing something wrong?

If the starting value for b is in the interval 0.4, 0.6, then the
fitted curve looks good. I guess a rough starting value for the
frequency is necessary.

My guess is that the curvature of the objective function is only
*locally* convex in the frequency, in b.
Also, it looks to me there is an indeterminacy in c (modulo 2*pi, or
it looks like modulo pi with switch in sign for b)

That's the problems that leastsq/curvefit might have, estimating in
frequency domain might be more appropriate. I guess somebody has a
better solution.

Josef

>
> any ideas, or references to places I can read about it, would be great!
>
>                thanks,
>
>                        Brian Blais
>
> --
> Brian Blais
> bblais@bryant.edu
> http://web.bryant.edu/~bblais
> http://bblais.blogspot.com/
>
> #=====================code below====================
> import pylab
> from scipy import optimize
> from numpy import *
>
> def m1(t,a,b,c):
>    return a*cos(b*t+c)
>
>
> x=pylab.linspace(0,50,100)
> y=cos(x/2)+pylab.randn(len(x))*.1
>
> pylab.figure(1)
> pylab.clf()
> pylab.plot(x,y,'o-')
>
> func=m1
> p0=[5,0.9,2]
>
> popt, pcov = optimize.curve_fit(func, x, y,
>        p0=p0,  # initial guess
>        )
>
>
> xf=linspace(x[0],x[-1],1000)
> yf=func(xf,*popt)
>
> pylab.plot(xf,yf,'.-')
>
> pylab.draw()
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list