为什么 scipy.optimize.curve_fit 不适合阻尼正弦曲线？

Question

注意：这是这个问题的后续

用户 Nick ODell 对于 curve_fit 不能很好地拟合正态曲线的问题给出了一个很好的答案——事实证明，相对于频率的目标函数并不“好”，并且有很多局部最小值与全局最小值相差甚远。为了解决这个问题，您应该首先使用 FFT 近似频率，并将其作为您的初始猜测。这是有效的，因为它将优化器置于全局最优附近，因此它不会在其他地方“丢失”。

现在，我有另一个数据集，这次是阻尼正弦曲线。我尝试使用相同的 FFT 技巧来获得良好的曲线拟合，但我得到的只是一条直线。

我的代码在这个pastebin中给出（如你所见，我发现

optimize.basinhopping

工作得很好，但我希望有更好的适合）。

没有长长的数据列表的代码：

ydata=[float(_) for _ in ydata.strip().splitlines()]
 
 
from scipy.optimize import curve_fit as cf
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import basinhopping as bh
xdata=np.arange(0.05,500.01,0.05)
 
def damped_sin_fun(x,a,b,c,d,e):
    return a*(np.exp(-x/e))*np.sin(b*x+c)+d
 
fft = np.fft.fft(ydata - np.array(ydata).mean())
timestep = xdata[1] - xdata[0]
freq = np.fft.fftfreq(len(fft), d=timestep)
largest_component = np.abs(fft).argmax()
phase_guess = np.angle(fft[largest_component]) * freq[largest_component]
initial_frequency_guess = freq[largest_component] * 2 * np.pi
 
p_opt,p_cov=cf(damped_sin_fun,xdata,ydata, p0=[0.05, initial_frequency_guess, np.pi/2 ,np.mean(ydata),2.3])
 
print(initial_frequency_guess)
print(p_opt)
#print(((damped_sin_fun(xdata,*p_opt)-ydata)**2).mean())
optimize=bh(lambda args: ((damped_sin_fun(xdata,*args)-ydata)**2).mean(), [0.05, initial_frequency_guess, np.pi/2 ,np.mean(ydata),2.3]).x
 
 
plt.plot(xdata,damped_sin_fun(xdata,*p_opt))
plt.xlim((0,10))
#plt.plot(xdata,damped_sin_fun(xdata,*optimize))
#plt.plot(xdata,ydata, 'r.-', ms=1)
plt.show()

Answer 1

我会通过将 Y 分解为三个信号来解决这个问题：指数项、常数项和正弦项。

y = exp * sine + constant

然后，我会分别安装这三个，然后再将它们组合起来。

为了拟合指数项，您不能只将其拟合到信号，因为您还不知道正弦参数应该是什么。然而，您可以做的是选择一个指数项，使标准差沿着信号保持恒定。一旦信号沿其长度具有恒定幅度，就可以使用上一个问题中描述的过程来拟合正弦波。

但是，人们无法从字面上找到每一点的标准差。相反，我所做的是计算滚动窗口中的标准差，窗口大小为 100。（这是以样本为单位，而不是时间单位。）对于具有 100 个样本进行平均之前的点，它使用向后充满。这是一个有点任意的参数，但我发现 20 到 1000 之间的值效果相当好。

samples_per_window = 100
exp_trend_data = pd.Series(ydata).rolling(window=samples_per_window).std().bfill().values

拟合指数项的完整代码：

import pandas as pd
import numpy as np
from scipy.optimize import curve_fit
from matplotlib import pyplot as plt


xdata = np.arange(0.05,500.01,0.05)
ydata = np.loadtxt('test579_data.txt')

samples_per_window = 100
exp_trend_data = pd.Series(ydata).rolling(window=samples_per_window).std().bfill().values

def exp_term(x, amplitude1, rate):
    return amplitude1 * np.exp(-x / rate)

beginning_value = exp_trend_data[:100].mean()
end_value = exp_trend_data[-100:].mean()
decrease_ratio = beginning_value / end_value
p0 = [beginning_value, np.max(xdata) / np.log(decrease_ratio)]
p_opt_exp, _ = curve_fit(exp_term, xdata, exp_trend_data, p0=p0)

plt.plot(exp_term(xdata, *p_opt_exp))
plt.plot(exp_trend_data)
plt.yscale('log')

这是指数项（蓝色）与目标（橙色）的图表。

（注意：这里的 Y 轴是对数的。在半对数图上，任何指数项都是直线。但是，目标不是直线。这意味着这个问题中的衰减不是严格的指数。）

接下来，拟合常数项。

constant_term = ydata.mean()

接下来，拟合正弦项。通过重新排列方程

y = exp * sine + constant

，我们可以单独得到正弦项。

sine_term_data = (ydata - constant_term) / exp_term(xdata, *p_opt_exp)
plt.plot(sine_term_data)

剧情：

这并没有完全消除标准差的变化，但已经足够接近了。

接下来，拟合正弦项。

def find_dominant_frequency(xdata, ydata):
    """Find dominiant frequency in units of radians"""
    fft = np.fft.fft(ydata - np.array(ydata).mean())
    timestep = xdata[1] - xdata[0]
    freq = np.fft.fftfreq(len(fft), d=timestep)
    fft = fft[:len(fft) // 2]
    freq = freq[:len(freq) // 2]
    largest_component = np.abs(fft).argmax()
    phase_guess = np.angle(fft[largest_component]) * freq[largest_component]
    initial_frequency_guess = freq[largest_component] * 2 * np.pi
    return initial_frequency_guess, phase_guess


def sine_term(x, amplitude2, freq, phase):
    return amplitude2*np.sin(freq*x+phase)


initial_frequency_guess, phase_guess = find_dominant_frequency(xdata, sine_term_data)
p0 = [np.sqrt(2), initial_frequency_guess, phase_guess]
p_opt_sine, _ = curve_fit(sine_term, xdata, sine_term_data, p0=p0)
plt.plot(xdata, sine_term_data)
plt.plot(xdata, sine_term(xdata, *p_opt_sine))
pos = 300
plt.xlim(pos, pos+20)

接下来，所有术语可以组合在一起。

def damped_sin_fun(x, amplitude, freq, phase, intercept, rate):
    return amplitude * np.exp(-x / rate) * np.sin(freq * x + phase) + intercept

amplitude1, rate = p_opt_exp
amplitude2, freq_guess, phase_guess = p_opt_sine
p0 = [
    amplitude1 * amplitude2,
    freq_guess,
    phase_guess,
    constant_term,
    rate,
]
plt.plot(xdata, damped_sin_fun(xdata, *p0))
plt.plot(xdata, ydata)

您可以运行另一条曲线拟合来一次调整所有参数。这与独立解决方案的结果非常相似，但略有不同。

p_opt, _ = curve_fit(damped_sin_fun, xdata, ydata, p0=p0)
print(p0)
print(p_opt.tolist())
plt.plot(xdata, ydata)
plt.plot(xdata, damped_sin_fun(xdata, *p_opt))
pos = 300
plt.xlim(pos, pos+20)

这会让你感到不适。