如何使用 pylab 和 numpy 将正弦曲线拟合到我的数据中？

Posted 2023-03-12

【中文标题】如何使用 pylab 和 numpy 将正弦曲线拟合到我的数据中？

【英文标题】：How do I fit a sine curve to my data with pylab and numpy?

【发布时间】：2013-05-18 22:52:28

【问题描述】：

我试图表明经济体遵循相对正弦的增长模式。我正在构建一个 python 模拟，以表明即使我们让某种程度的随机性占据主导地位，我们仍然可以产生相对正弦的东西。

我对我生成的数据感到满意，但现在我想找到一些方法来获得与数据非常匹配的正弦图。我知道您可以进行多项式拟合，但您可以进行正弦拟合吗？

【参考方案1】：

以上所有答案都基于曲线拟合，并且大多数使用迭代方法 - 它们都工作得很好，但我想使用 FFT 添加不同的方法。在这里，我们转换数据，将除峰值频率之外的所有数据都设置为零，然后进行逆变换。请注意，您可能希望在执行 FFT 之前删除数据均值（和去趋势），然后您可以在之后将它们添加回来。

import numpy as np
import pylab as plt

# fake data

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
f = 1.05
data = 3.0*np.sin(f*t+0.001) + np.random.randn(N) # create artificial data with noise

# FFT...
mfft=np.fft.fft(data)
imax=np.argmax(np.absolute(mfft))
mask=np.zeros_like(mfft)
mask[[imax]]=1
mfft*=mask
fdata=np.fft.ifft(mfft)


plt.plot(t, data, '.')
plt.plot(t, fdata,'.', label='FFT')
plt.legend()
plt.show()

【参考方案2】：

您可以在 scipy 中使用 least-square optimization 函数将任意函数拟合到另一个函数。在拟合 sin 函数的情况下，要拟合的 3 个参数是偏移 ('a')、幅度 ('b') 和相位 ('c')。

只要您对参数提供合理的初步猜测，优化应该会很好地收敛。幸运的是，对于正弦函数，其中 2 个的初步估计很容易：可以通过取数据的平均值和通过 RMS (3*标准偏差/sqrt(2)) 获得的幅度。

注意：作为稍后的编辑，还添加了频率拟合。这不能很好地工作（可能导致非常糟糕的配合）。因此，请自行决定使用，我的建议是不要使用频率拟合，除非频率误差小于几个百分点。

这导致以下代码：

import numpy as np
from scipy.optimize import leastsq
import pylab as plt

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
f = 1.15247 # Optional!! Advised not to use
data = 3.0*np.sin(f*t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise

guess_mean = np.mean(data)
guess_std = 3*np.std(data)/(2**0.5)/(2**0.5)
guess_phase = 0
guess_freq = 1
guess_amp = 1

# we'll use this to plot our first estimate. This might already be good enough for you
data_first_guess = guess_std*np.sin(t+guess_phase) + guess_mean

# Define the function to optimize, in this case, we want to minimize the difference
# between the actual data and our "guessed" parameters
optimize_func = lambda x: x[0]*np.sin(x[1]*t+x[2]) + x[3] - data
est_amp, est_freq, est_phase, est_mean = leastsq(optimize_func, [guess_amp, guess_freq, guess_phase, guess_mean])[0]

# recreate the fitted curve using the optimized parameters
data_fit = est_amp*np.sin(est_freq*t+est_phase) + est_mean

# recreate the fitted curve using the optimized parameters

fine_t = np.arange(0,max(t),0.1)
data_fit=est_amp*np.sin(est_freq*fine_t+est_phase)+est_mean

plt.plot(t, data, '.')
plt.plot(t, data_first_guess, label='first guess')
plt.plot(fine_t, data_fit, label='after fitting')
plt.legend()
plt.show()

编辑：我假设您知道正弦波中的周期数。如果您不这样做，则安装起来会有些棘手。您可以尝试通过手动绘图来猜测周期数，并尝试将其优化为您的第 6 个参数。

【讨论】：

这个解决方案虽然被 OP 接受，但似乎跳过了最棘手的部分：f 中的 频率 y = Amplitude*sin(frequency*x +Phase) + Offset。如果f 未知，此方法效果如何？

有理由不使用scipy的curve_fit函数吗？我猜最初的猜测和/或curve_fit对函数做出假设是有问题的。

我认为您提供函数的初始参数值的顺序是错误的。 est_a, est_b, est_c = leastsq(optimize_func, [guess_b, guess_a, guess_c])[0] 呢？为清楚起见，我建议将 _a 替换为 _offset，将 _b 替换为 _amp，将 _c 替换为 _phase，并在 lambda 中使用 x[i] 的递增顺序。

【参考方案3】：

这是一个不需要手动猜测频率的无参数拟合函数fit_sin()：

import numpy, scipy.optimize

def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return "amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)

初始频率猜测由使用 FFT 的频域中的峰值频率给出。假设只有一个主频（零频峰值除外），拟合结果几乎是完美的。

import pylab as plt

N, amp, omega, phase, offset, noise = 500, 1., 2., .5, 4., 3
#N, amp, omega, phase, offset, noise = 50, 1., .4, .5, 4., .2
#N, amp, omega, phase, offset, noise = 200, 1., 20, .5, 4., 1
tt = numpy.linspace(0, 10, N)
tt2 = numpy.linspace(0, 10, 10*N)
yy = amp*numpy.sin(omega*tt + phase) + offset
yynoise = yy + noise*(numpy.random.random(len(tt))-0.5)

res = fit_sin(tt, yynoise)
print( "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )

plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt, yynoise, "ok", label="y with noise")
plt.plot(tt2, res["fitfunc"](tt2), "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

即使噪音很大，结果也很好：

幅度=1.00660540618，角频率=2.03370472482，相位=0.360276844224，偏移=3.95747467506，最大。 Cov.=0.0122923578658

【讨论】：

亲爱的 unsym 我试图运行您的代码，但不幸的是我收到以下消息：当我尝试绘制函数时，TypeError: 'numpy.float64' 对象不能被解释为整数。你有什么办法解决这个问题吗？

我刚刚按原样运行代码，它使用 python 3.6 运行没有任何错误（在笔记本中打开 matplotlib“内联”绘图的 jupyter 笔记本中）

【参考方案4】：

对我们更友好的是函数curvefit。举个例子：

import numpy as np
from scipy.optimize import curve_fit
import pylab as plt

N = 1000 # number of data points
t = np.linspace(0, 4*np.pi, N)
data = 3.0*np.sin(t+0.001) + 0.5 + np.random.randn(N) # create artificial data with noise

guess_freq = 1
guess_amplitude = 3*np.std(data)/(2**0.5)
guess_phase = 0
guess_offset = np.mean(data)

p0=[guess_freq, guess_amplitude,
    guess_phase, guess_offset]

# create the function we want to fit
def my_sin(x, freq, amplitude, phase, offset):
    return np.sin(x * freq + phase) * amplitude + offset

# now do the fit
fit = curve_fit(my_sin, t, data, p0=p0)

# we'll use this to plot our first estimate. This might already be good enough for you
data_first_guess = my_sin(t, *p0)

# recreate the fitted curve using the optimized parameters
data_fit = my_sin(t, *fit[0])

plt.plot(data, '.')
plt.plot(data_fit, label='after fitting')
plt.plot(data_first_guess, label='first guess')
plt.legend()
plt.show()

【讨论】：

@IceArdor：你能添加一个你建议的求解器的工作代码示例吗？

@Vasco 如果我有一个 m*n 形状的 x_train 和 m*1 形状的 y_train 那么在这种情况下我会收到这个错误：ValueError：操作数不能与形状一起广播（38563 ,54) (38563,) 。那我现在该怎么办？

@Vasco 你能给我看一些例子吗，或者如果可以在这里用一些示例代码发布你的答案，那将非常有帮助：***.com/questions/57027040/…

【参考方案5】：

当前将 sin 曲线拟合到给定数据集的方法需要首先猜测参数，然后是一个交互过程。这是一个非线性回归问题。

由于方便的积分方程，另一种方法是将非线性回归转换为线性回归。这样就不需要初始猜测，也不需要迭代过程：直接得到拟合。

如果是函数y = a + r*sin(w*x+phi) 或y=a+b*sin(w*x)+c*cos(w*x)，请参阅"Régression sinusoidale" 发表在Scribd 上的论文"Régression sinusoidale" 的第35-36 页

在函数y = a + p*x + r*sin(w*x+phi) 的情况下：“混合线性和正弦回归”一章的第 49-51 页。

对于更复杂的函数，一般流程在章节"Generalized sinusoidal regression"第54-61页进行了说明，后面是数值示例y = r*sin(w*x+phi)+(b/x)+c*ln(x)，第62-63页

【讨论】：

碰巧，你有没有关于阻尼正弦波变换的参考资料？

+1 这看起来很有趣，问题是参考文献是法语的。有谁知道一些英语资源（可能是博客文章/论文/代码中的任何内容），可以更详细地说明该方法？如果确实可以转换为线性回归，我想知道为什么这种方法不受欢迎？

啊，我找到了一个很棒的资源：@​​987654322@。它不仅提供了实现，还附带了整篇论文的translation。