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

【问题标题】：Error in nonlinear least square fit of sine function with exponential decay具有指数衰减的正弦函数的非线性最小二乘拟合误差
【发布时间】：2020-10-22 06:23:17
【问题描述】：

我的非线性数据使用公式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, 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()

【问题讨论】：

请阅读如何创建 MCVE：stackoverflow.com/help/minimal-reproducible-example 我们喜欢在 Stackoverflow 上提供帮助，但我们大多数人没有太多时间提供帮助。借助 MCVE，您可以让我们轻松重现您的问题，以便我们可以花更多时间专注于您遇到的实际问题。
请提供准确的输入和准确的预期输出。也可以链接到理论。
查看您的数据：与幅度相比，您的衰减很大，导致小值问题。您有大约 10 个数据点的行为类似于您的预期公式，但有 90 个点基本上为零 - 最重要的是，您用随机噪声掩盖了这些值。您的方法适用于较小的衰减值和较大的安培值。

标签： python lmfit

【解决方案1】：

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

如果你绘制了你在拟合中实际使用的数据：

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

你会看到这个：

尽管您实际上使用了differential_evolution，但您还将获得的最差最佳匹配标记为leastsq。

如果你降低了测试数据中的噪声和衰减，并且真正符合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, yn, 'ko', 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

还有的情节

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

【讨论】：