Error in nonlinear least square fit of sine function with exponential decay

Question

My nonlinear data are approximated using a least square fit with the formula Asin(wt+phase)exp(-decay*t) while keeping omega(w) as a constant. I have tried several approaches without success.

Below is my code

import numpy as np
from numpy import loadtxt
import lmfit
np.random.seed(2)
x = np.linspace(0, 10, 101)
decay = 4.5
shift = 0
amp = 0.0015
y = amp * np.sin(x*5+shift) * np.exp(-x*decay)
yn = y + np.random.normal(size=y.size, scale=0.450)


def resid(params, x, ydata):
    decay = params['decay'].value
    shift = params['shift'].value
    amp = params['amp'].value

    y_model = amp * np.sin(x*5+shift) * np.exp(-x*decay)
    return y_model - ydata
params = lmfit.Parameters()
params.add('shift', 0.0, min=-np.pi, max=np.pi)
params.add('amp', 0.0015, min=0, max=0.02)
params.add('decay', 4.0, min=0, max=10.0)
fit = lmfit.minimize(resid, params, args=(x, yn), method='differential_evolution')
print("\n\n# Fit using differential_evolution:")
lmfit.report_fit(fit)
plt.plot(x, y, 'ko', lw=2)
plt.plot(x, yn+fit.residual, 'b--', lw=2)
plt.legend(['data', 'leastsq', 'diffev'], loc='upper left')
plt.show()

Answer 1

Comments pointed out that the decay you have is very severe. But also, the noise you add to yn is enormous compared to your decaying sine wave. You and others may have missed this because you do not plot the yn array that you fit, rather you plot y, without the added noise data.

If you do plot the data you actually use in the fit:

plt.plot(x, yn, 'ko', lw=2)
plt.show()

you will see this:

You also label the very poort best-fit you get as leastsq though you actually used differential_evolution.

If you lower the noise and the decay in your test data, and really fit with leastsq, you might have something like this:

import numpy as np
import lmfit
import matplotlib.pyplot as plt

np.random.seed(2)
x = np.linspace(0, 10, 101)

decay = 0.6
shift = 0
amp = 0.025
y = amp * np.sin(x*5+shift) * np.exp(-x*decay)
yn = y + np.random.normal(size=y.size, scale=0.001)

def resid(params, x, ydata):
    decay = params['decay'].value
    shift = params['shift'].value
    amp = params['amp'].value

    y_model = amp * np.sin(x*5+shift) * np.exp(-x*decay)
    return y_model - ydata

params = lmfit.Parameters()
params.add('shift', 0.0, min=-np.pi, max=np.pi)
params.add('amp', 0.01)
params.add('decay', 0.2, min=0, max=50)
fit = lmfit.minimize(resid, params, args=(x, yn))


print("\n\n# Fit using differential_evolution:")
lmfit.report_fit(fit)
plt.plot(x, yn, 'ko', lw=2)
plt.plot(x, yn+fit.residual, 'b--', lw=2)
plt.legend(['data', 'leastsq'], loc='upper right')
plt.show()

which will give a report of

# Fit using differential_evolution:
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 21
    # data points      = 101
    # variables        = 3
    chi-square         = 1.0894e-04
    reduced chi-square = 1.1117e-06
    Akaike info crit   = -1381.71974
    Bayesian info crit = -1373.87437
[[Variables]]
    shift:  0.00796328 +/- 0.01999031 (251.03%) (init = 0)
    amp:    0.02448871 +/- 7.5756e-04 (3.09%) (init = 0.01)
    decay:  0.59826922 +/- 0.02587107 (4.32%) (init = 0.2)
[[Correlations]] (unreported correlations are < 0.100)
    C(amp, decay)   =  0.725
    C(shift, amp)   = -0.147
    C(shift, decay) = -0.105

and a plot of

I think the conclusion is that your residual function and setting up the problem are OK, but your very noisy test data set is not represented by that decaying sine-wave very well.