具有指数衰减的正弦函数的非线性最小二乘拟合中的误差

问题描述

我的非线性数据使用最小二乘拟合与公式Asin(wt+phase)exp(-decay*t)进行近似，同时将omega(w)保持不变。我尝试了几种方法都没有成功。

下面是我的代码

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,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()

解决方法

评论指出，您的衰减非常严重。而且，与衰减的正弦波相比，您添加到yn的噪声非常大。您和其他人可能错过了这一点，因为您没有绘制自己适合的yn数组，而是绘制了y，而没有添加噪声data。

如果您绘制了实际用于拟合的数据：

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

您将看到以下内容：

尽管您实际上使用的是leastsq，但您也将最差的最佳拟合标记为differential_evolution。

如果您降低测试数据中的噪声和衰减，并且确实适合leastsq，则可能会出现以下情况：

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,lw=2)
plt.plot(x,yn+fit.residual,'b--',lw=2)
plt.legend(['data','leastsq'],loc='upper right')
plt.show()

将给出报告

# 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

和情节

我认为结论是您的残差函数和设置问题都可以，但是您的非常嘈杂的测试数据集不能很好地由衰减正弦波表示。