Parameters estimation for curve fitting with Python lmfit

Question

I am using Python lmfit to do least square fit with monthly average data from 2005-2016. I have constructed the function like below: equation and the original code shown as below:

# t is in fractional years, e.g. 2017+122./365.
def fun(t, a, b, c, A1, A2, A3, A4, B1, B2, B3, B4):
    An=[A1,A2,A3,A4]
    Bn=[B1,B2,B3,B4]
    sum=np.sum([An[i] * np.sin(2 * np.pi * (i + 1) * t+Bn[i]) for i in range(len(An))])
    return a+b*t+c*t*t+sum

mod = Model(fun)
pars = mod.make_params(a=-10, b=0.003, c=0.01, A1=-1., A2=1., A3=1., A4=1., B1=-1., B2=1., B3=1., B4=1.)

result = mod.fit(y, pars, t=t)
print(result.fit_report())

plt.plot(t, y, 'bo')
plt.plot(t, result.best_fit, 'r-')
plt.show()

fitted line and the original data dots It seems that the Fourier terms didn't work. Therefore, I am curious that how to give a suitable initial estimation on the function parameters such as A1, A2,A3...?

Moreover, you should provide minimal working example, i.e. at least a complete working code including all your imports. A minimal data set would be nice as well. — mikuszefski
Apart from that, a can be checked from mean. The Bs are probably best set to zero. A1 to something like 0.5 * ( max() - min() ). The other A's to zero. (Assuming that all your data looks similar to what you show in the link), b and c to zero as well. — mikuszefski
You might try using the scipy module scipy.optimize.differental_evolution to determine the initial parameters, I have an example of using this scipy module to fit a fouble Lorentzian peak equation to Raman spectroscopy of carbon nanotubes at bitbucket.org/zunzuncode/ramanspectroscopyfit — James Phillips
Thanks for helpful comments from mikuszefski and James. I made mistakes on the shape of the sum over the Fourier terms. See answer by M Newville. And see the fit plot at lh3.googleusercontent.com/-qE0-_ojL5N4/Wi85OooEJPI/AAAAAAAAGq8/… — Cheng

M Newville M Newville · Accepted Answer · 2017-12-11T23:22:55

np.sum does not do what you want it to do. It will sum your expression to a single scalar value, not an array of the same length as t. That scalar value then collapses your parameters A1, ... B4 onto a single value, and the fit will have no way to determine these values.

I think you want to make a 2D array of shape (4, len(t)) and then sum over only the first dimension, leaving an array of len(t) that is the sum over the 4 Fourier component.

Try replacing your

sum=np.sum([An[i]*np.sin(2*np.pi*(i+1)*t+Bn[i]) for i in range(len(An))])

with

sum=np.array([An[i]*np.sin(2*np.pi*(i+1)*t+Bn[i]) for i in range(len(An))]).sum(axis=0)

Parameters estimation for curve fitting with Python lmfit

1 Answers