在将正弦曲线拟合到周期数据时，如何改进scipy.Optimize.curve_fit参数的初始猜测，否则如何改进拟合？
问题描述
我有一个包含测量的space-separated csv文件。第一列是测量时间，第二列是相应的测量值，第三列是误差。The file can be found here.我想使用Python将函数g的参数a_i, f, phi_n与数据进行匹配：
正在读取数据：
import numpy as np
data=np.genfromtxt('signal.data')
time=data[:,0]
signal=data[:,1]
signalerror=data[:,2]
绘制数据：
import matplotlib.pyplot as plt
plt.figure()
plt.plot(time,signal)
plt.scatter(time,signal,s=5)
plt.show()
获取结果：
现在让我们计算一下周期信号的初步频率猜测：
from gatspy.periodic import LombScargleFast
dmag=0.000005
nyquist_factor=40
model = LombScargleFast().fit(time, signal, dmag)
periods, power = model.periodogram_auto(nyquist_factor)
model.optimizer.period_range=(0.2, 10)
period = model.best_period
我们得到结果：0.5467448186001437
我将该函数定义为适合N=10：
def G(x, A_0,
A_1, phi_1,
A_2, phi_2,
A_3, phi_3,
A_4, phi_4,
A_5, phi_5,
A_6, phi_6,
A_7, phi_7,
A_8, phi_8,
A_9, phi_9,
A_10, phi_10,
freq):
return (A_0 + A_1 * np.sin(2 * np.pi * 1 * freq * x + phi_1) +
A_2 * np.sin(2 * np.pi * 2 * freq * x + phi_2) +
A_3 * np.sin(2 * np.pi * 3 * freq * x + phi_3) +
A_4 * np.sin(2 * np.pi * 4 * freq * x + phi_4) +
A_5 * np.sin(2 * np.pi * 5 * freq * x + phi_5) +
A_6 * np.sin(2 * np.pi * 6 * freq * x + phi_6) +
A_7 * np.sin(2 * np.pi * 7 * freq * x + phi_7) +
A_8 * np.sin(2 * np.pi * 8 * freq * x + phi_8) +
A_9 * np.sin(2 * np.pi * 9 * freq * x + phi_9) +
A_10 * np.sin(2 * np.pi * 10 * freq * x + phi_10))
现在我们需要一个适合G的函数：
def fitter(time, signal, signalerror, LSPfreq):
from scipy import optimize
pfit, pcov = optimize.curve_fit(lambda x, _A_0,
_A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit:
G(x, _A_0, _A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit),
time, signal, p0=[11, 2, 0, #p0 is the initial guess for numerical fitting
1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq],
sigma=signalerror, absolute_sigma=True)
error = [] # DEFINE LIST TO CALC ERROR
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i]) ** 0.5) # CALCULATE SQUARE ROOT OF TRACE OF COVARIANCE MATRIX
except:
error.append(0.00)
perr_curvefit = np.array(error)
return pfit, perr_curvefit
检查我们获得的内容：
LSPfreq=1/period
pfit, perr_curvefit = fitter(time, signal, signalerror, LSPfreq)
plt.figure()
model=G(time,pfit[0],pfit[1],pfit[2],pfit[3],pfit[4],pfit[5],pfit[6],pfit[7],pfit[8],pfit[8],pfit[10],pfit[11],pfit[12],pfit[13],pfit[14],pfit[15],pfit[16],pfit[17],pfit[18],pfit[19],pfit[20],pfit[21])
plt.scatter(time,model,marker='+')
plt.plot(time,signal,c='r')
plt.show()
收益率：
这显然是错误的。如果我玩函数fitter定义中的初始猜测p0，我可以得到稍微好一点的结果。设置
p0=[11, 1, 0,
0.1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq]
给我们(放大)：
哪一个好一点。高频分量仍然存在，尽管高频分量的幅度被猜测为零。根据对数据的目测，原始的p0似乎也比修改后的版本更合理。
我尝试了p0的不同值，虽然更改p0确实会更改结果，但我得不到与数据很好匹配的线条。
为什么此模型拟合方法失败？我如何才能改进Get Good Fit？
The whole code can be found here.
此问题的原始版本已发布here。
编辑：
PyCharm对代码的p0部分发出警告：
预期类型为‘Union[None，int，Float，Complex]’，实际类型为‘List[Union[int，Float]，Any]]’
我不知道如何处理，但可能与此相关。
解决方案
对于在噪声数据上计算最佳拟合周期模型，典型的基于优化的方法通常在除最人为的情况外的所有情况下都会失败。这是因为成本函数在频率空间将是高度多峰的，因此任何没有密集网格搜索的优化方法几乎肯定会陷入局部最小。
在这种情况下，最佳密集网格搜索将是您用来查找初始值的Lomb-Scarger周期图的变体，您可以跳过优化步骤，因为Lomb-Scargle已经为您优化了它。
目前最好的广义Lomb-ScarglePython实现可在Astropy中获得(完全披露：此实现的大部分都是我编写的)。上面使用的模型在那里称为截断傅立叶模型，可以通过为nterms参数指定适当的值来进行拟合。
使用您的数据，您可以从拟合和绘制具有五个傅立叶项的广义周期图开始：
from astropy.stats import LombScargle
ls = LombScargle(time, signal, signalerror, nterms=5)
freq, power = ls.autopower()
plt.plot(freq, power);
freq, power = ls.autopower(maximum_frequency=24)
plt.plot(freq, power);
这向我们展示的是网格上每个频率下最适合的傅立叶模型的有效逆卡方。 我们现在可以找到最佳频率并计算该频率的最佳模型：
best_freq = freq[power.argmax()]
tfit = np.linspace(time.min(), time.max(), 10000)
signalfit = ls.model(tfit, best_freq)
plt.errorbar(time, signal, signalerror, fmt='.k', ecolor='lightgray');
plt.plot(tfit, signalfit)
plt.xlim(time[500], time[800]);
如果您对模型参数本身感兴趣，则可以使用Lomb-Scarger算法背后的低级例程。
from astropy.stats.lombscargle.implementations.mle import design_matrix
X = design_matrix(time, best_freq, signalerror, nterms=5)
parameters = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, signal / signalerror))
print(parameters)
# [ 1.18351382e+01 2.24194359e-01 5.72266632e-02 -1.23807286e-01
# 1.25825666e-02 7.81944277e-02 -1.10571718e-02 -5.49132878e-02
# 9.51544241e-03 3.70385961e-02 9.36161528e-06]
这些是线性化模型的参数，即
signal = p_0 + sum_{n=1}^{5}[p_{2n - 1} sin(2pi n f t) + p_{2n} cos(2pi n f t)]
这些线性sin/cos幅值可以通过一点三角运算转换回非线性幅值和相位。
我相信这将是将多项傅里叶级数与您的模型相匹配的最佳方法，因为它避免了对表现不佳的成本函数进行优化，并使用快速算法使基于网格的计算变得容易处理。
相关文章