scipy.optimize on high frequency sine function

Question

I am using Python 2.7. I am wondering why the optimize function of SciPy doesn't converge to the right function when the target is a high frequency sinus wave.

import numpy as np
from scipy import optimize

test_func = lambda x: 5*np.sin(15*x+3)+1
t = linspace(0,25,100000)
y_t = test_func(t)
plot(t,y_t)

fitfunc = lambda p, x: p[0]*np.sin(p[1]*x+p[2])+p[3]
errfunc = lambda p, x, y: fitfunc(p, x) - y 

p0 = [max(y_t),10,2,0]
p1, success = optimize.leastsq(errfunc, p0, args=(t,y_t))


plot(t,fitfunc(p1,t))

One can clearly see that the end solution diverges from the target clearly. Am I doing something wrong ? Is the error function ill adapted here ?

Thanks for any input

Answer 1

Your problem is that there are a large number of local minima in your residuals function as the phase and frequency shift with respect to their true values; without really good initial guesses for the phase and frequency you will converge into one instead of falling into the much deeper, global minimum:

If you don't have any more information about the phase and frequency, you can either estimate them from a FFT of the data or rewrite your formula as

Asin(bx + phi) + d = Acos(phi)sin(bx) + Asin(phi)cos(bx) + d

which has only one nonlinear parameter (b): you can use a grid-search for b and much faster and more reliable linear least-squares fitting for the rest (a1 = Acos(phi), a2 = Asin(phi) and d).

Here's a plot of the rms residual as the frequency, b varies, showing the various minima: