2

I am trying to fit a simple sinewave to data and applying bounds to help constrain the fit.

I am using the 2nd answer posted here which is working perfectly, however when I apply bounds as per convention, e.g.:

def fit_sin(tt, yy):
    import scipy.optimize
    import numpy as np
    '''
    Fit sin to the input time sequence, and return dict of fitting parameters:
    "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"
    '''
    tt = np.array(tt)
    yy = np.array(yy)
    ff = np.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(np.fft.fft(yy))
    guess_freq = abs(ff[np.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = np.std(yy) * 2.**0.5
    guess_offset = np.mean(yy)
    guess = np.array([guess_amp, 2.*np.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * np.sin(w*t + p) + c
    boundary = ([-np.inf, -np.inf, -np.pi, -np.inf],[np.inf, np.inf, np.pi, np.inf])
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess, bounds=boundary)
    A, w, p, c = popt
    f = w/(2.*np.pi)
    fitfunc = lambda t: A * np.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": np.max(pcov), "rawres": (guess,popt,pcov)}

I am getting the following error:

RuntimeError: Optimal parameters not found: The maximum number of function evaluations is exceeded.

If I set the bounds to infinity and I change the boundary line to simply,

boundary = ([-np.inf, -np.inf, -np.inf, -np.inf],[np.inf, np.inf, np.inf, np.inf])

the function works.

To test the function you can do:

import pylab as plt

N, amp, omega, phase, offset, noise = 500, 1., 2., .5, 4., 3
#N, amp, omega, phase, offset, noise = 50, 1., .4, .5, 4., .2
#N, amp, omega, phase, offset, noise = 200, 1., 20, .5, 4., 1
tt = numpy.linspace(0, 10, N)
tt2 = numpy.linspace(0, 10, 10*N)
yy = amp*numpy.sin(omega*tt + phase) + offset
yynoise = yy + noise*(numpy.random.random(len(tt))-0.5)

res = fit_sin(tt, yynoise)
print( "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )

plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt, yynoise, "ok", label="y with noise")
plt.plot(tt2, res["fitfunc"](tt2), "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

In particular I am just trying to constrain the phase p of the fit to be between -np.pi and np.pi, the others can be infinity.

The actual data I'm using this function on is clean and the fit is quick and really accurate, but it's sometimes fitting the phase a whole cycle 'out of step' when the data starts at or near +/- np.pi, it will fit it 2*np.pi out phase.

I can't manually catch this either as it's fitting to many thousands of data sets and I'm looking at relative phase differences between them all.

Any help appreciated, thanks.

rh1990
  • 880
  • 7
  • 17
  • 32
  • What is the `boundary` line that is giving you troubles? If I use this: `boundary = ([0.0, 0.0, -np.pi, -100.0],[10.0, 1000.0, np.pi, 100.0])` the code works as expected. Note that the fitting algorithm is different depending on whether you provide bounds or not. As a side note, `import ...` statements are typically outside of the function, the docstring is enclosed with triple `"""` double-quotes, and immediately follows the function definition. – norok2 Oct 08 '18 at 14:29
  • 1
    Thanks for your input. I've just found the issue: the fit was actually fitting to negative values of amplitude e.g. `-100` instead of `100`, and then that mean the fitted phase value was forced to be `pi` larger. – rh1990 Oct 08 '18 at 14:47
  • The code in the question right now works well. Hence, unless you also include the code that was causing you troubles, I am unsure it would be of any value to future readers. – norok2 Oct 08 '18 at 15:03

1 Answers1

1

I fixed it, for a particular data set causing me trouble, my function was fitting negative amplitudes which then forced my phase value to be np.pi larger than anticipated.

I'll leave this answer here for anyone having similar troubles with sine-wave fitting!

rh1990
  • 880
  • 7
  • 17
  • 32