我有一个包含测量space-separated csv
文件。第一列是测量的时间,第二列是相应的测量值,第三列是错误。 The file can be found here.我想以适应函数g的参数a_i, f, phi_n
到数据,使用Python:
在读取数据:
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()
产量:
这显然是错误的。如果我在功能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]
让我们(放大):
这是一个好一点。 ,高频成分仍然是尽管高频compontents的幅度被猜测是零。原来p0
似乎也比基于数据的目视检查修改后的版本更加合理。
我有不同的价值观起到周围p0
,虽然确实在改变p0
改变的结果,我没有得到一个合理的线很好拟合数据。
为什么这个模型拟合方法失败?我怎样才能提高获得一个不错的选择?
The whole code can be found here.
这个问题的原始版本被张贴here。
编辑:
PyCharm给出了代码的p0
部分的警告:
预计式“联盟[无,整数,浮点,复杂]”,得到了“一览[联盟[整型,浮点],任何]”,而不是
与我不知道该如何处理,但可能是相关的。
在有噪声的数据计算最适合的周期性模式,典型的基于优化的方法通常会失败,但在所有的最做作的情况。这是因为成本函数就会有强烈的多模态频率空间,所以任何优化方法短期密集网格搜索的几乎肯定会陷入局部最小。
在这种情况下,最好的密集网格搜索将是你用来寻找初始值博士伦-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);
混叠作用是立即清除这里:由于您的数据点之间的间隔,上述约24所有频率都只是在考虑频率低于24下使用该信号的别名,让我们重新计算只是周期图的相关部分:
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]);
如果你有兴趣在模型参数本身,有可能获得在使用低级别的例程博士伦-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)]
这些线性正弦/余弦的振幅可以被转换回非线性振幅和相位与位三角的。
我相信这将是拟合多项傅里叶级数到模型的最好的方法,因为它避免了在优化不良的行为成本的功能,并采用快速算法,使基于网格的计算驯服。
这里是出现给一个确定的数据拟合的代码。这使用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()