如何改进参数的初始猜测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()

Answer 2

对于计算噪声数据下的最佳拟合周期模型，典型的基于优化的方法通常在所有情况下都会失败，但大多数人为的情况除外。这是因为代价函数在频率空间中是高度多模态的，因此任何一种缺少密集网格搜索的优化方法几乎肯定会陷入局部极小。在

在这种情况下，最佳密集网格搜索将是用于查找初始值的Lomb Scargle周期图的变体，您可以跳过优化步骤，因为Lomb Scargle已经为您优化了它。在

通用lombscargle的最好的Python实现可以在Astropy中找到（完全公开：我编写了这个实现的大部分内容）。上面使用的模型称为截断Fourier模型，可以通过为nterms参数指定适当的值来进行拟合。在

使用您的数据，您可以从拟合和绘制带有五个傅立叶项的广义周期图开始：

from astropy.stats import LombScargle
ls = LombScargle(time, signal, signalerror, nterms=5)
freq, power = ls.autopower()
plt.plot(freq, power);

这里的混叠很明显：由于数据点之间的间隔，所有高于24的频率都只是频率低于24的信号的别名。考虑到这一点，让我们重新计算周期图的相关部分：

^{pr2}$

这向我们展示了在网格上每一个频率下，最佳拟合傅里叶模型的卡方倒数。 我们现在可以找到最佳频率并计算该频率的最佳拟合模型：

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

如果您对模型参数本身感兴趣，可以使用lombscargle算法后面的低级例程。在

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(2\pi n f t) + p_{2n} cos(2\pi n f t)]

这些线性正弦/余弦振幅可以转换回非线性振幅和相位与一点三角学。在

我相信这将是你的模型拟合一个多项傅立叶级数的最佳方法，因为它避免了对性能较差的成本函数进行优化，并使用快速算法使基于网格的计算更容易处理。在