哪些技术可以有效地找到任意数据点的周期性？
What techniques are effective to find periodicity in arbitrary data points?
我所说的“任意”是指我没有在适合进行 FFT 的网格上采样的信号。我只有事件发生的时间点（例如时间），我想要估计的速率，例如：
p = [0, 1.1, 1.9, 3, 3.9, 6.1 ...]
...可能是标称周期（重复间隔）为 1.0 的过程命中，但有噪声和一些漏检。
是否有处理此类数据的众所周知的方法？
听起来你需要决定你到底想确定什么。如果你想知道一组时间戳中的平均间隔，那很容易（取平均值或中位数）。
如果您预计间隔可能会发生变化，那么您需要了解它的变化速度。然后你可以找到一个 windowed 移动平均线。您需要知道它变化的速度有多快，以便您可以 select 您的 window 大小适当 - 较大的 window 会给您更平滑的结果，但较小的 [=22] =] 将对 faster-changing 比率更加敏感。
如果您不知道数据是否遵循某种模式，那么您可能正处于数据探索领域。在那种情况下，我将从绘制间隔开始，看看是否出现了一种模式。如果数据非常嘈杂，这也可能受益于应用移动平均线。
基本上，数据中是否包含某些内容及其含义取决于您和您对该领域的了解。也就是说，在任何一组时间戳中， 都会 是平均值（您也可以轻松计算方差以指示数据的可变性），但这取决于您是否这个平均值没有任何意义。
如果正确初始化，最小二乘算法可以解决问题。为此可以应用聚类方法。
执行 FFT 时，信号被描述为正弦波的总和。频率的幅度可以描述为由least square fit on the signal产生的。因此，如果信号被不均匀采样，如果要估计傅立叶变换，解决相同的最小二乘问题可能是有意义的。如果应用于均匀采样的信号，它归结为相同的结果。
由于您的信号是离散的，您可能希望将其作为 Dirac combs. It seems more sound to minimize the sum of squared distance to the nearest Dirac of the Dirac comb. This is a non-linear optimization problem where Dirac combs are described by their period and offset. This non-linear least-square problem can be solved by mean of the Levenberg-Marquardt algorithm. Here is an python example making use of the scipy.optimize.leastsq() function. Moreover, the error on the estimated period and offset can be estimated as depicted in How to compute standard deviation errors with scipy.optimize.least_squares . It is also documented in the documentation of curve_fit() and Getting standard errors on fitted parameters using the optimize.leastsq method in python
的总和来拟合
然而，周期的一半，或周期的三分之一，......，也适合，并且周期的倍数是局部最小值，可以通过优化 [=17= 的初始化来避免] 已应用。
请注意，该过程可以扩展到多维数据以寻找具有不同基本周期的周期模式或混合周期模式。
import numpy as np
from scipy.optimize import least_squares
from scipy.optimize import leastsq
from sklearn.cluster import MeanShift, estimate_bandwidth
ticks=[0,1.1,1.9,3,3.9,6.1]
import scipy
print scipy.__version__
def crudeEstimate():
# loooking for the period by looking at differences between values :
diffs=np.zeros(((len(ticks))*(len(ticks)-1))/2)
k=0
for i in range(len(ticks)):
for j in range(i):
diffs[k]=ticks[i]-ticks[j]
k=k+1
#see
X = np.array(zip(diffs,np.zeros(len(diffs))), dtype=np.float)
bandwidth = estimate_bandwidth(X, quantile=1.0/len(ticks))
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X)
labels = ms.labels_
cluster_centers = ms.cluster_centers_
print cluster_centers
labels_unique = np.unique(labels)
n_clusters_ = len(labels_unique)
for k in range(n_clusters_):
my_members = labels == k
print "cluster {0}: {1}".format(k, X[my_members, 0])
estimated_period=np.min(cluster_centers[:,0])
return estimated_period
def disttoDiracComb(x):
residual=np.zeros((len(ticks)))
for i in range(len(ticks)):
mindist=np.inf
for j in range(len(x)/2):
offset=x[2*j+1]
period=x[2*j]
#print period, offset
index=np.floor((ticks[i]-offset)/period)
#print 'index', index
currdist=ticks[i]-(index*period+offset)
if currdist>0.5*period:
currdist=period-currdist
index=index+1
#print 'event at ',ticks[i], 'not far from index ',index, '(', currdist, ')'
#currdist=currdist*currdist
#print currdist
if currdist<mindist:
mindist=currdist
residual[i]=mindist
#residual=residual-period*period
#print x, residual
return residual
estimated_period=crudeEstimate()
print 'crude estimate by clustering :',estimated_period
xp=np.array([estimated_period,0.0])
#res_1 = least_squares(disttoDiracComb, xp,method='lm',xtol=1e-15,verbose=1)
p,pcov,infodict,mesg,ier=leastsq(disttoDiracComb, x0=xp,ftol=1e-18, full_output=True)
#print ' p is ',p, 'covariance is ', pcov
# see
s_sq = (disttoDiracComb(p)**2).sum()/(len(ticks)-len(p))
pcov=pcov *s_sq
perr = np.sqrt(np.diag(pcov))
#print 'estimated standard deviation on parameter :' , perr
print 'estimated period is ', p[0],' +/- ', 1.96*perr[0]
print 'estimated offset is ', p[1],' +/- ', 1.96*perr[1]
应用于您的示例，它打印：
crude estimate by clustering : 0.975
estimated period is 1.0042857141346768 +/- 0.04035792507868619
estimated offset is -0.011428571139828817 +/- 0.13385206912205957
我所说的“任意”是指我没有在适合进行 FFT 的网格上采样的信号。我只有事件发生的时间点（例如时间），我想要估计的速率，例如：
p = [0, 1.1, 1.9, 3, 3.9, 6.1 ...]
...可能是标称周期（重复间隔）为 1.0 的过程命中，但有噪声和一些漏检。
是否有处理此类数据的众所周知的方法？
听起来你需要决定你到底想确定什么。如果你想知道一组时间戳中的平均间隔，那很容易（取平均值或中位数）。
如果您预计间隔可能会发生变化，那么您需要了解它的变化速度。然后你可以找到一个 windowed 移动平均线。您需要知道它变化的速度有多快，以便您可以 select 您的 window 大小适当 - 较大的 window 会给您更平滑的结果，但较小的 [=22] =] 将对 faster-changing 比率更加敏感。
如果您不知道数据是否遵循某种模式，那么您可能正处于数据探索领域。在那种情况下，我将从绘制间隔开始，看看是否出现了一种模式。如果数据非常嘈杂，这也可能受益于应用移动平均线。
基本上，数据中是否包含某些内容及其含义取决于您和您对该领域的了解。也就是说，在任何一组时间戳中， 都会 是平均值（您也可以轻松计算方差以指示数据的可变性），但这取决于您是否这个平均值没有任何意义。
如果正确初始化，最小二乘算法可以解决问题。为此可以应用聚类方法。
执行 FFT 时，信号被描述为正弦波的总和。频率的幅度可以描述为由least square fit on the signal产生的。因此，如果信号被不均匀采样，如果要估计傅立叶变换，解决相同的最小二乘问题可能是有意义的。如果应用于均匀采样的信号，它归结为相同的结果。
由于您的信号是离散的，您可能希望将其作为 Dirac combs. It seems more sound to minimize the sum of squared distance to the nearest Dirac of the Dirac comb. This is a non-linear optimization problem where Dirac combs are described by their period and offset. This non-linear least-square problem can be solved by mean of the Levenberg-Marquardt algorithm. Here is an python example making use of the scipy.optimize.leastsq() function. Moreover, the error on the estimated period and offset can be estimated as depicted in How to compute standard deviation errors with scipy.optimize.least_squares . It is also documented in the documentation of curve_fit() and Getting standard errors on fitted parameters using the optimize.leastsq method in python
然而，周期的一半，或周期的三分之一，......，也适合，并且周期的倍数是局部最小值，可以通过优化 [=17= 的初始化来避免] 已应用。
请注意，该过程可以扩展到多维数据以寻找具有不同基本周期的周期模式或混合周期模式。
import numpy as np
from scipy.optimize import least_squares
from scipy.optimize import leastsq
from sklearn.cluster import MeanShift, estimate_bandwidth
ticks=[0,1.1,1.9,3,3.9,6.1]
import scipy
print scipy.__version__
def crudeEstimate():
# loooking for the period by looking at differences between values :
diffs=np.zeros(((len(ticks))*(len(ticks)-1))/2)
k=0
for i in range(len(ticks)):
for j in range(i):
diffs[k]=ticks[i]-ticks[j]
k=k+1
#see
X = np.array(zip(diffs,np.zeros(len(diffs))), dtype=np.float)
bandwidth = estimate_bandwidth(X, quantile=1.0/len(ticks))
ms = MeanShift(bandwidth=bandwidth, bin_seeding=True)
ms.fit(X)
labels = ms.labels_
cluster_centers = ms.cluster_centers_
print cluster_centers
labels_unique = np.unique(labels)
n_clusters_ = len(labels_unique)
for k in range(n_clusters_):
my_members = labels == k
print "cluster {0}: {1}".format(k, X[my_members, 0])
estimated_period=np.min(cluster_centers[:,0])
return estimated_period
def disttoDiracComb(x):
residual=np.zeros((len(ticks)))
for i in range(len(ticks)):
mindist=np.inf
for j in range(len(x)/2):
offset=x[2*j+1]
period=x[2*j]
#print period, offset
index=np.floor((ticks[i]-offset)/period)
#print 'index', index
currdist=ticks[i]-(index*period+offset)
if currdist>0.5*period:
currdist=period-currdist
index=index+1
#print 'event at ',ticks[i], 'not far from index ',index, '(', currdist, ')'
#currdist=currdist*currdist
#print currdist
if currdist<mindist:
mindist=currdist
residual[i]=mindist
#residual=residual-period*period
#print x, residual
return residual
estimated_period=crudeEstimate()
print 'crude estimate by clustering :',estimated_period
xp=np.array([estimated_period,0.0])
#res_1 = least_squares(disttoDiracComb, xp,method='lm',xtol=1e-15,verbose=1)
p,pcov,infodict,mesg,ier=leastsq(disttoDiracComb, x0=xp,ftol=1e-18, full_output=True)
#print ' p is ',p, 'covariance is ', pcov
# see
s_sq = (disttoDiracComb(p)**2).sum()/(len(ticks)-len(p))
pcov=pcov *s_sq
perr = np.sqrt(np.diag(pcov))
#print 'estimated standard deviation on parameter :' , perr
print 'estimated period is ', p[0],' +/- ', 1.96*perr[0]
print 'estimated offset is ', p[1],' +/- ', 1.96*perr[1]
应用于您的示例，它打印：
crude estimate by clustering : 0.975
estimated period is 1.0042857141346768 +/- 0.04035792507868619
estimated offset is -0.011428571139828817 +/- 0.13385206912205957