[SciPy-User] fitting cosine with curve_fit - problem with frequency
josef.pktd@gmai...
josef.pktd@gmai...
Wed Oct 12 20:09:59 CDT 2011
On Wed, Oct 12, 2011 at 8:55 PM, <josef.pktd@gmail.com> wrote:
> 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.
some guessing, but it works in the examples I tried:
import matplotlib.pyplot as plt
psd, fr = plt.psd(y)
midx = np.argmax(psd)
print fr[midx], psd[midx]
p0[1] = fr[midx]*np.pi*2
Josef
>
> 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