尝试使用 scipy 将三角函数拟合到数据

Question

我正在尝试使用 scipy.optimize.curve_fit 拟合一些数据。我有 read the documentation and also this Whosebug post，但似乎都没有回答我的问题。

我有 some data，它是简单的二维数据，看起来很像三角函数。我想用一般的三角函数来适应它使用 scipy.

我的做法如下：

from __future__ import division
import numpy as np
from scipy.optimize import curve_fit



#Load the data
data = np.loadtxt('example_data.txt')
t = data[:,0]
y = data[:,1]


#define the function to fit
def func_cos(t,A,omega,dphi,C):
    # A is the amplitude, omega the frequency, dphi and C the horizontal/vertical shifts
    return A*np.cos(omega*t + dphi) + C

#do a scipy fit
popt, pcov = curve_fit(func_cos, t,y)

#Plot fit data and original data
fig = plt.figure(figsize=(14,10))
ax1 = plt.subplot2grid((1,1), (0,0))

ax1.plot(t,y)
ax1.plot(t,func_cos(t,*popt))

这输出：

其中蓝色是数据，橙色是拟合。显然我做错了什么。有什么指点吗？

Answer 1

如果没有为参数 p0 的初始猜测提供值，则假定每个参数的值为 1。来自文档：

p0 : array_like, optional
Initial guess for the parameters (length N). If None, then the initial values will all be 1 (if the number of parameters for the function can be determined using introspection, otherwise a ValueError is raised).

由于您的数据具有非常大的 x 值和非常小的 y 值，1 的初始猜测与实际解决方案相去甚远，因此优化器不会收敛。您可以通过提供可以从数据中猜测/近似的合适的初始参数值来帮助优化器：

振幅：A = (y.max() - y.min()) / 2
偏移量：C = (y.max() + y.min()) / 2
频率：这里我们可以通过将连续的y值相乘来估计过零的数量，并检查哪些乘积小于零。这个数字除以总 x 范围给出了频率，为了以 pi 为单位得到它，我们可以将该数字乘以 pi：y_shifted = y - offset; oemga = np.pi * np.sum(y_shifted[:-1] * y_shifted[1:] < 0) / (t.max() - t.min())
相移：可设为零，dphi = 0

所以综上所述，可以使用以下初始参数猜测：

offset = (y.max() + y.min()) / 2
y_shifted = y - offset
p0 = (
    (y.max() - y.min()) / 2,
    np.pi * np.sum(y_shifted[:-1] * y_shifted[1:] < 0) / (t.max() - t.min()),
    0,
    offset
)
popt, pcov = curve_fit(func_cos, t, y, p0=p0)

这给了我以下拟合函数：

尝试使用 scipy 将三角函数拟合到数据

Trying to fit a trig function to data with scipy

python

curve-fitting

scipy

scipy-optimize