【发布时间】：2019-06-26 11:45:31
【问题描述】：
我有一个包含测量值的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 似乎也比基于数据视觉检查的修改版本更合理。
我对@987654348@ 使用了不同的值，虽然更改p0 确实会改变结果，但我没有得到一条与数据相当吻合的线。
为什么这种模型拟合方法会失败？我怎样才能更好地适应？
The whole code can be found here.
这个问题的原始版本发布在here。
编辑：
PyCharm 对代码的p0 部分给出警告：
预期类型 'Union[None,int,float,complex]'，得到了 'List[Union[int,float],Any]]'
我不知道如何处理，但可能是相关的。
【问题讨论】：
-
我有一个建议。首先尝试使用较小的数据子集，例如“time=time[:500];signal=signal[:500];signalerror=signalerror[:500]”，因为它运行得更快、更容易跟...共事。拟合数据子集的结果可以作为整个数据集的初始参数估计。
标签： python scipy curve-fitting