如何改进参数的初始猜测scipy.optimize.curve U拟合在将正弦拟合到周期性数据时，还是改进拟合？

Question

我有一个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]

绘制数据：

^{pr2}$

得到结果：

现在让我们计算周期信号的初步频率猜测：

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确实会改变结果，但我没有得到一条与数据相当吻合的行。在

为什么这种模型拟合方法失败？我怎样才能更好的适应？

The whole code can be found here.

---

这个问题的原始版本已经发布了here。在

---

编辑：

PyCharm对代码的p0部分发出警告：

Expected type 'Union[None,int,float,complex]', got 'List[Union[int,float],Any]]' instead

我不知道怎么处理，但可能有关系。在

Answer 1

下面的代码看起来与数据匹配良好。这使用scipy的差分进化（DE）遗传算法来估计曲线拟合（）的初始参数。为了加速遗传算法，代码使用前500个数据点的数据子集进行初始参数估计。虽然结果看起来很好，但是这个问题有一个复杂的错误空间，有很多参数，而且遗传算法需要一些时间来运行（在我的史前笔记本上大约15分钟）。您应该考虑在午餐时间或晚上使用完整的数据集进行测试，以验证拟合的参数是否有任何有用的改进。DE的scipy实现使用拉丁超立方体算法来确保对参数空间的彻底搜索，这需要搜索范围-请检查示例的边界是否合理。在

import numpy as np

from scipy.optimize import differential_evolution
import warnings

data=np.genfromtxt('signal.data')

time=data[:,0]
signal=data[:,1]
signalerror=data[:,2]

# value for reduced size data set used in initial parameter estimation
# to sllow the genetic algorithm to run faster than with all data
geneticAlgorithmSlice = 500

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
LSPfreq=1/period


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))



# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = G(time[:geneticAlgorithmSlice], *parameterTuple)
    return np.sum((signal[:geneticAlgorithmSlice] - val) ** 2.0)

def generate_Initial_Parameters():
    parameterBounds = []
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([-50.0, 50.0])
    parameterBounds.append([LSPfreq/2.0, LSPfreq*2.0])

    # "seed" the numpy random number generator for repeatable results
    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x


print("Starting genetic algorithm...")
# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()
print("Genetic algorithm completed")


def fitter(time, signal, signalerror, initialParameters):

    from scipy import optimize

    pfit, pcov = optimize.curve_fit(G, time, signal, p0=initialParameters,
                                    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


pfit, perr_curvefit = fitter(time, signal, signalerror, geneticParameters)

plt.figure()
model=G(time,*pfit) 
plt.scatter(time,model,marker='+')
plt.plot(time,model)
plt.plot(time,signal,c='r')
plt.show()