我怎样才能提高参数的初始猜测scipy.optimize.curve_fit装修时正弦周期性数据,或以其他方式提高配件?

问题描述 投票:2回答:2

我有一个包含测量space-separated csv文件。第一列是测量的时间,第二列是相应的测量值,第三列是错误。 The file can be found here.我想以适应函数g的参数a_i, f, phi_n到数据,使用Python:

enter image description here

在读取数据:

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

获得的结果:

enter image description here

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

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

产量:

enter image description here

这显然是错误的。如果我在功能p0定义初始猜测fitter玩,我可以得到一个稍微好一点的结果。设置

p0=[11,  1, 0,
        0.1, 0,
        0, 0,
        0, 0,
        0, 0,
        0, 0,
        0, 0,
        0, 0,
        0, 0,
        0, 0,
        LSPfreq]

让我们(放大):

enter image description here

这是一个好一点。 ,高频成分仍然是尽管高频compontents的幅度被猜测是零。原来p0似乎也比基于数据的目视检查修改后的版本更加合理。

我有不同的价值观起到周围p0,虽然确实在改变p0改变的结果,我没有得到一个合理的线很好拟合数据。

为什么这个模型拟合方法失败?我怎样才能提高获得一个不错的选择?

The whole code can be found here.


这个问题的原始版本被张贴here


编辑:

PyCharm给出了代码的p0部分的警告:

enter image description here

预计式“联盟[无,整数,浮点,复杂]”,得到了“一览[联盟[整型,浮点],任何]”,而不是

与我不知道该如何处理,但可能是相关的。

python scipy curve-fitting
2个回答
3
投票

在有噪声的数据计算最适合的周期性模式,典型的基于优化的方法通常会失败,但在所有的最做作的情况。这是因为成本函数就会有强烈的多模态频率空间,所以任何优化方法短期密集网格搜索的几乎肯定会陷入局部最小。

在这种情况下,最好的密集网格搜索将是你用来寻找初始值博士伦-Scargle周期图的一个变种,你可以跳过优化步骤,因为博士伦-Scargle已经优化它。

广义博士伦-Scargle最好的Python的实现目前在Astropy缴费(全面披露:我写的散装该实现的)。你以上使用该模型被称为那里作为截断的傅里叶模型,并且可以通过指定nterms参数的适当的值来适应。

使用你的数据,你可以通过安装并开始密谋五傅立叶项广义周期图:

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

enter image description here

混叠作用是立即清除这里:由于您的数据点之间的间隔,上述约24所有频率都只是在考虑频率低于24下使用该信号的别名,让我们重新计算只是周期图的相关部分:

freq, power = ls.autopower(maximum_frequency=24)
plt.plot(freq, power);

enter image description here

什么这表明我们实际上是在网格上的每个频率的最佳拟合模型傅立叶逆卡方。现在,我们可以找到最佳的频率并计算该频率的最适合的模式:

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

enter image description here

如果你有兴趣在模型参数本身,有可能获得在使用低级别的例程博士伦-scargle算法落后。

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

这些线性正弦/余弦的振幅可以被转换回非线性振幅和相位与位三角的。

我相信这将是拟合多项傅里叶级数到模型的最好的方法,因为它避免了在优化不良的行为成本的功能,并采用快速算法,使基于网格的计算驯服。


1
投票

这里是出现给一个确定的数据拟合的代码。这使用SciPy的的差分演化(DE)的遗传算法来估计初始参数curve_fit()。为了加快遗传算法,该代码使用用于初始参数估计第一500个数据点的数据子集。虽然结果在视觉上看起来不错,这个问题有许多参数的复杂的误差空间,遗传算法将需要一些时间来运行(在我的笔记本电脑史前近15分钟)。你应该使用完整数据的午餐时间或深夜设置验证拟合参数是否没有任何有用的改进考虑测试。该SciPy的实现DE的使用拉丁超立方体算法,以确保彻底的搜索参数空间,这就需要在其中边界搜索 - 请检查该示例的边界似乎是合理的。

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()
© www.soinside.com 2019 - 2024. All rights reserved.