为您推荐:
精华内容
最热下载
问答
  • 5星
    17KB younow22 2021-02-15 20:33:07
  • 5星
    1.61MB m0_54656292 2021-08-26 23:28:23
  • 5星
    14KB weixin_42696271 2021-09-11 03:41:54
  • 5星
    43KB qq_34093397 2021-01-04 10:45:38
  • 5星
    11.8MB guoruibin123 2021-05-21 15:45:32
  • 5星
    10KB weixin_42696333 2021-09-10 22:42:23
  • 5星
    5KB m0_53867522 2021-05-25 12:53:29
  • 5星
    145.04MB qq_31988139 2021-08-10 09:54:59
  • 5星
    23.83MB HowardEmily 2021-01-14 11:46:42
  • 5星
    296KB m0_52957036 2020-09-06 23:56:19
  • 1.25MB weixin_38569203 2021-04-05 11:48:09
  • 837KB weixin_38577200 2021-02-26 03:31:08
  • 2KB weixin_38629042 2021-05-28 20:35:52
  • 针对特定的预测问题,只是拥有数据还不够,想要从纷繁复杂的数据关系中挖掘出可用于预测的规律或模式,还得运用恰当的分析方法。比如聚类分析,恰当地选择聚类算法,可以按维度将数据适当地分群,根据各类的特征制订...

            针对特定的预测问题,只是拥有数据还不够,想要从纷繁复杂的数据关系中挖掘出可用于预测的规律或模式,还得运用恰当的分析方法。比如聚类分析,恰当地选择聚类算法,可以按维度将数据适当地分群,根据各类的特征制订营销计划或决策,抑或是根据各类不冋规律建立起更有针对性的预测模型;还有常用的关联分析,可以从事物的历史数据中挖掘出变化规律有指导性地对未来进行预测,如此等等。本内容将分别介绍常用的分析方法,本节介绍相关分析。

    目录

    1. 自相关分析            2. 偏相关分析           3. 简单相关分析           4. 互相关分析           5. 典型相关分析


           相关关系是一种与函数关系相区别的非确定性关系,而相关分析就是研究事物或现象之间是否存在这种非确定性关系的统计方法。相关分析按处理问题的不同,通常可分为自相关分析、偏相关分析、简单相关分析、互相关分析及典型相关分析。其中:

    • 自相关分析、偏相关分析:适用于分析变量自身的规律;
    • 简单相关分析:通常可分析任意两个等长数列间的相关性;
    • 互相关分析:允许在一定的间隔下进行简单相关分析;
    • 典型相关分析:适用于分析两组变量的相关性。

    变量与变量的关系通常有两种:

    • 函数关系:表示确定的非随机变量之间的关系,可以用表达式具体地定义。比如,位移与速度就是确定的函数关系。
    • 非确定性关系:表示非随机变量与随机变量的关系,也就是说当给定一个变量时,另一个变量是随机的,它不能由具体的表达式来定义。比如,身高和体重就无法用一个确定的函数关系来表达,但是根据经验,身高比较高的人,其体重也比较重,说明两者是相关关系

    1. 自相关分析

            自相关是指同一时间序列在不同时刻取值的相关程度,假设有时间序列X_{t}, t=1,2,3,...,则在时刻t和滞后n阶t+n之间的相关即为n阶自相关,其定义如下: 

                                

    其中,函数f为计算相关系数的函数,可通过上式计算滞后n阶自相关系数的值。这里使用 airmiles时序数据来分析时间序列的自相关性,该数据集记录的是从1937年到1960年美国商业航空公司的飞机里程数据,如图1-1所示。 

                           

                                                      图 1-1 从1937年到1960年美国商业航空公司的飞机里程趋势

    在R语言中,可用acf函数分析时间序列的自相关性,该函数定义及参数说明如下表: 

    acf(airmiles,type='correlation',lag.max=10)
    
    # airmiles:美国1937-1960年客运里程营收(实际售出机位乘以飞行哩数)

    效果图如下: 

                       

    如图可见,滞后阶数为0时,相关系数为1,随着滞后阶数的增加,相关系数逐渐减弱,并趋于稳定。 

    2. 偏相关分析 

           在分析当期与前n期指标的相关性时,1期到n-1期之间的指标会对相关性的分析有影响,于是将这些影响指标去掉后,再进行的相关分析叫作偏相关分析。

            由于自相关性分析的是时间序列X_{t}, t=1,2,3,...,在时刻tt+n之间取值的相关性程度,其值是未在限定X_{t+1}\sim X_{t+n-1}取值的情况下进行计算的,所得的自相关系数多少会受X_{t+1}\sim X_{t+n-1}取值的影响。为了更加真实地计算自相关关系数值,需要在限定其他值的情况下进行计算,这就是所谓的偏相关,其定义如下:

            

          其中pf函数是求解X_{t}X_{t+n}在排除X_{t+1}\sim X_{t+n-1}因素影响的情况下的偏相关系数。同时,对X_{s}X_{t}在限定X_{k}的情况下,其偏相关系数定义如下: 

                                           

         在R语言中,可直接使用pacf函数分析时间序列的偏相关性,该函数定义及参数说明如表3-1-2所示 :

             

    pacf(airmiles,lag.max=10)

     效果图如下: 

     

           如图所示,最小为1阶滞后,对应值为0.876,与对应的1阶自相关系数相等,随着滞后阶数的增加(大于2阶),偏相关系数一直较小并且稳定。

    1阶滞后:相当于X_{t}X_{t+1},因此其自相关系数和偏自相关系数相等。

    3. 简单相关分析

           相关关系是一种非确定的关系,就好像身高与体重的关系一样,它们之间不能用一个固定的函关系来表示。而相关分析就是研究这种随机变量间相关关系的统计方法。此处,主要探讨不同特征对研究对象的相关性影响。常见进行相关分析的方法,主要有散点图和相关图。

    3.1 散点图

           散点图就是数据点在直角坐标系上的分布图,通常分为散点图矩阵和三维散点图。其中散点矩阵是由变量两两组合由数据点分布图构成的矩阵,而三维散点图就是从所有变量中选择三个变量进行绘制,进一步在三维空间里观察数据的形态。

    (1)散点图矩阵

           常用于绘制散点矩阵的方法是R语言中:

    • car包中的 scatterplotMatrix函数;
    • graphics包中的pais函数。

           这里以R语言自带的iris数据集为例,说明分析鸢尾花的四个指标的相关关系,绘制散点图矩阵,代码如下。

    # 图3-1
    pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris, 
       main="Simple Scatterplot Matrix")
    
    # 图3-2
    library(car)
    scatterplotMatrix(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width|Species, 
           data=iris)
    
    # 图3-3
    library(rgl)
    scatter3d(iris$Sepal.Length, iris$Petal.Length, iris$Petal.Width)

                

                                                                                          图 3-1

           如图3-1所示为所有变量的两两组合的散点图矩阵,每个散点图中呈现的是任意两变量的数据点,可通过数据点的分布,了解变量之间的相关性。此图中 Petal.Length与 Petal.Width对应的散点图比较接近线性,说明这两个变量间的相关性较强。

               

                                                                                                     图 3-2

           如代码所示,通过“|”指定了分组的变量,这里使用鸢尾花的种类进行分组。如图3-2所示对角线上的图形表示各个变量在不同鸢尾花类型下的分布情况;其他图形分别用不同颜色为数据点着色,同时实线表示对数据点的拟合线,上下虚线表示浮动范围,直接表示线性拟合线。根据该图可以更进一步地知道不同类型鸢尾花各变量的相关关系,以及线性及非线性的变化规律

    (2)三维散点图

            常用于绘制三维散点图的方法是:

    • car包的 scatter3d函数

            通过该函数绘制的图形,可以进行转动,以便切换不同角度,查看数据在三维空间的分布情况。通过使用 scatter3d函数绘制鸢尾花的 Sepal.Length、 Petal.Length、 Petal.Width这三个指标在三维空间的散点图,如下所示,该函数为三维空间中的点拟合了线性平面,通过旋转可以更直观地观察数据的分布规律。

                                             

                                                                                                图 3-3

    3.2 相关图

            所谓相关图是基于变量间的相关关系大小,通过可视化方式反应不同变量组合间相关关系的差异的图形。可以把相关图分为相关矩阵图、相关层次图。

    (1)相关矩阵图

           R语言中,绘制相关矩阵图的包,主要有两个,分别为:

    • corrgram包中的 corrgram函数;

            panel:统一设置非对角的面板,

            lower.panel和 upper panel:分别设置相关矩阵下三角和上三角的面板形状,

            通常可以设置:panel.pie(饼图)、 panel.shade(阴影图)、 panel.ellipse(椭圆图)、panel.pts(散点图)

    library(corrgram)
    
    #1、设置排序处理
    corrgram(mtcars,order=TRUE)
    
    #2、设置上下三角面板形状
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie)
    
    #3、只显示下三角部分
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=NULL)
    
    #4、调整面板颜色
    corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,
             col.regions=colorRampPalette(c("darkgoldenrod4","burlywood1","white",
             "darkkhaki","darkgreen")))

     

                                                                    图 3-4:以上程序1、2、3、4效果图

    • complot包中的 complot函数

             method:可以取7个值,分别为 circle(圆形,默认) 、square(方形) 、ellipse(椭圆)、number(数字)、pie(饼图) 、shade(阴影) 、color(热力图)。

            type:用来设置相关矩阵图的显示区域,通常取full全部,默认)upper(上三角)、lower(下三角)

            col:可通过 colorRampPalette函数向col参数赋值来设置颜色

            order:用来指定相关系数变量的排序方式,通常可取 original(原始顺序,默认),AOE(特征向量的角顺序,FPC(第一主成分的顺序), hclust(层次聚类的顺序)以及 alphabet(字母表顺序)

    library(corrplot)
    
    #1、设置method=color绘制热力矩阵图
    corrplot(cor(mtcars), method="color", order = "AOE",tl.col="black",tl.srt=45,
             addCoef.col="black",col=colorRampPalette(c("#7F0000","red","#FF7F00",
             "yellow","white", "cyan", "#007FFF", "blue","#00007F"))(20))
    
    #2、绘制上下三角及不同色彩的相关矩阵图
    library(RColorBrewer)
    par(mfrow=c(2,2))
    corrplot(cor(mtcars),type="lower")
    corrplot(cor(mtcars),type="lower",order="hclust", col=brewer.pal(n=8,name="RdYlBu"))
    corrplot(cor(mtcars),type="upper",order="AOE", col=c("black","white"),bg="lightblue")
    corrplot(cor(mtcars),type="upper",order="FPC", col=brewer.pal(n=8, name="PuOr"))
    par(mfrow=c(1,1))

        

                                                                         图 3-5:以上程序1、2效果图 

    (2)相关层次图

          层次图,通过计算样本间的距离来判断各样本是否属于同一类。通过将相关系数转化为距离度量,进行系统聚类,旨在分析各变量的相关关系及组合影响情况。通常有四种方法将相关系数转化为相异性度量d

     使用mtcars数据集,进行系统聚类,如下: 方式一:hclust() 

    • d=1-r
    • d=(1-r)/2
    • d=1-|r|
    • d=\sqrt{1-r^{2}}
    d<-sqrt(1-cor(mtcars)^2)
    hc<-hclust(as.dist(d))
    plot(hc)
    rect.hclust(hc,k=3)  # 添加聚类分类矩形,如分为3类

     效果如下图所示: 

                           

    方式二:pvclust(),通过p值来评估层次聚类的不确定性,表明了影响显著的聚类,显著的类越多,效果越好。

    参数:method.dist() 指定用于计算距离的方法,设定为correlation时,按 d=1-r 计算;设定为abscor时,按 d=1-|r| 计算。

    library(pvclust)
    cluster.bootstrap <- pvclust(mtcars, nboot=1000, method.dist="correlation")
    plot(cluster.bootstrap)
    pvrect(cluster.bootstrap)  # 添加聚类分类矩形

     效果如下图所示:

                         

    • AU值:Approxiamtely Unbiased,通过多尺度自助重抽样法计算而来;
    • BP值:Bootstrap Probability,通过普通自助重抽样法计算而来 

    4. 互相关分析

          与自相关不同,互相关是指两个时间序列在做任意两个不同时刻的相关程度,假设有时间序列X_{t}, t=1,2,3,...Y_{t}, t=1,2,3,...,则X在时刻t和y在时刻t+n之间的相关即为n阶互相关,其定义如下:

                              

            其中函数∫为计算相关系数的函数,可通过上式计算滞后n阶互相关系数的值。

            使用 airmiles和 Lakehuron数据集,来说明互相关分析的方法。其中 airmiles数据集记录了从1937年到1960年美国商业航空公司每年的乘客里程数据,而 Lakehuron数据集记录了从1875年到1972年休伦湖每年的湖平面的测量数据。将它们的时间限制在1937年到1960年,并绘制各自的时间序列曲线如图4-1所示。 

                                     

                          

                                                                                          图 4-1 

    在R语言中,可直接使用ccf函数分析时间序列的互相关性,该函数定义及参数说明如表3-1-3所示:

            

            

    plot(airmiles)
    plot(LakeHuron,xlim=c(1937,1960))
    
    ccf(airmiles, ts(LakeHuron, start = 1937, end = 1960), type="correlation")

    结果如下: 

                         

            效果如图所示,当没有延迟(即Lag=0)时,互相关系数取得最小,将近-0.8,两个时间序列大致上呈负相关关系。而当延迟10阶时,两个时间序列互相关系数最大,且为弱正相关。可通过互相关性的分析,构建用于预测的合适指标。 

    5. 典型相关分析

            典型相关是指两组变量间的相关关系,当然,它不是指对两组变量进行两两组合的简单相关,而是反映两组变量作为两个整体之间的相关性。所以典型相关的问题就是如何构建综合指标使其相关性最大。常见的做法是:

            从两组变量中提取一对线性组合,使其相关性最大;然后,从剩余的且与前面不相关的线性组合中再提取一对,使其具有最大的相关性,这样依次进行下去。所提取的线性组合就是典型变量,两组典型变量对应的相关系数就是典型相关系数。

            R语言程序包 stats中的 cancor函数,可以直接求解两个数值矩阵的典型相关系数,它的定义及参数说明如表3-1-4所示。

                  

             现使用iris数据集,分别将第1、2列和3、4列转换成数值矩阵,使用 cancor函数求解典型相关系数,代码如下。

    x<-as.matrix(iris[,1:2])
    y<-as.matrix(iris[,3:4])
    cancor(x,y)

      结果如下: 

                                     

           由结果可知,典型相关系数有两个,分别是0.9409690和0.1239369,说明第一组典型变量的相关性很强,第二组典型变量的相关性很弱。这种情况,通常第一个典型相关系数用于分析。

           xcof和ycof分别表示矩阵x和y中典型变量与原变量对应的系数矩阵,可通过标准化的数据矩阵进行乘法运算得到典型变量。

     

    展开全文
    qq_27586341 2019-06-18 14:35:11
  • 1.31MB weixin_38571449 2021-05-11 15:34:11
  • 4.35MB weixin_38724106 2021-02-09 11:35:54
  • 21KB Lidx0710 2021-10-01 17:36:46
  • 150) for item in ax.get_xticks()]) plt.show() 如上图所示,是将时间序列分割成了 20 个等长的时间段,然后计算每个时间窗口的互相关。这给了我们更细粒度的视角来观察信号的相互作用。例如,在第一个窗口内(第一...

    作者:Andrew Chung

    公众号:WealthQuant

    链接:

    https://www.zhihu.com/question/23525783/answer/956912446

    已获得作者授权转载

    量化两个时间序列之间的相关性可以从很多方向着手, 下面说说我的总结仅供参考(Python). 基于你的信号类型,你对信号作出的假设,以及你想要从数据中寻找什么样的同步性数据的目标,来决定使用那种相关性测量.

    • Person相关

    • 时间滞后互相关(TLCC)以及加窗的 TLCC

    • 动态时间扭曲(DTW)

    • 瞬时相位同步


    1. 皮尔逊相关 —— 最简单也是最好的方法

    Person相关可以衡量两个连续信号如何随时间共同变化,并且可以以数字 -1(负相关)、0(不相关)和 1(完全相关)表示出它们之间的线性关系。它很直观,容易理解,也很好解释。但是当使用皮尔逊相关的时候,有两件事情需要注意,它们分别是:第一,异常数据可能会干扰相关评估的结果;第二,它假设数据都是同方差的,这样的话,数据方差在整个数据范围内都是同质的。通常情况下,相关性是全局同步性的快照测量法。所以,它不能提供关于两个信号间方向性的信息,例如,哪个信号是引导信号,哪个信号是跟随信号。

    如下的代码加载的就是样本数据(它和代码位于同一个文件夹下),并使用 Pandas 和 Scipy 计算皮尔逊相关,然后绘制出了中值滤波的数据。

    import pandas as pd
    import numpy as np
    import matplotlib.pyplot as plt
    import seaborn as sns
    import scipy.stats as stats
    
    
    df = pd.read_csv('synchrony_sample.csv')
    overall_pearson_r = df.corr().iloc[0,1]
    print(f"Pandas computed Pearson r: {overall_pearson_r}")
    # 输出:使用 Pandas 计算皮尔逊相关结果的 r 值:0.2058774513561943
    
    
    r, p = stats.pearsonr(df.dropna()['S1_Joy'], df.dropna()['S2_Joy'])
    print(f"Scipy computed Pearson r: {r} and p-value: {p}")
    # 输出:使用 Scipy 计算皮尔逊相关结果的 r 值:0.20587745135619354,以及 p-value:3.7902989479463397e-51
    
    
    # 计算滑动窗口同步性
    f,ax=plt.subplots(figsize=(7,3))
    df.rolling(window=30,center=True).median().plot(ax=ax)
    ax.set(xlabel='Time',ylabel='Pearson r')
    ax.set(title=f"Overall Pearson r = {np.round(overall_pearson_r,2)}")
    plt.show()
    

    再次重申,所有的皮尔逊 r 值都是用来衡量全局同步的,它将两个信号的关系精简到了一个值当中。尽管如此,使用皮尔逊相关也有办法观察每一刻的状态,即局部同步性。计算的方法之一就是测量信号局部的皮尔逊相关,然后在所有滑动窗口重复该过程,直到所有的信号都被窗口覆盖过。由于可以根据你想要重复的次数任意定义窗口的宽度,这个结果会因人而异。在下面的代码中,我们使用 120 帧作为窗口宽度(4 秒左右),然后在下图展示出我们绘制的每一刻的同步结果。

    # 设置窗口宽度,以计算滑动窗口同步性
    r_window_size = 120
    # 插入缺失值
    df_interpolated = df.interpolate()
    # 计算滑动窗口同步性
    rolling_r = df_interpolated['S1_Joy'].rolling(window=r_window_size, center=True).corr(df_interpolated['S2_Joy'])
    f,ax=plt.subplots(2,1,figsize=(14,6),sharex=True)
    df.rolling(window=30,center=True).median().plot(ax=ax[0])
    ax[0].set(xlabel='Frame',ylabel='Smiling Evidence')
    rolling_r.plot(ax=ax[1])
    ax[1].set(xlabel='Frame',ylabel='Pearson r')
    plt.suptitle("Smiling data and rolling window correlation")
    

    总的来说,皮尔逊相关是很好的入门学习教程,它提供了一个计算全局和局部同步性的很简单的方法。但是,它不能提供信号动态信息,例如哪个信号先出现,而这个可以用互相关来衡量。


    2. 时间滞后互相关 —— 评估信号动态性

    时间滞后互相关(TLCC)可以定义两个信号之间的方向性,例如引导-追随关系,在这种关系中,引导信号会初始化一个响应,追随信号则重复它。还有一些其他方法可以探查这类关系,包括格兰杰因果,它常用于经济学,但是要注意这些仍然不一定能反映真正的因果关系。但是,通过查看互相关,我们还是可以提取出哪个信号首先出现的信息。

    如上图所示,TLCC 是通过逐步移动一个时间序列向量(红色线)并反复计算两个信号间的相关性而测量得到的。如果相关性的峰值位于中心(offset=0),那就意味着两个时间序列在此时相关性最高。但是,如果一个信号在引导另一个信号,相关性的峰值就可能位于不同的坐标值上。下面这段代码应用了一个使用了 pandas 提供功能的互相关函数。同时它也可以将数据打包,这样相关性边界值也能通过添加信号另一边的数据而计算出来。

    def crosscorr(datax, datay, lag=0, wrap=False):
        """ Lag-N cross correlation. 
        Shifted data filled with NaNs 
    
    
        Parameters
        ----------
        lag : int, default 0
        datax, datay : pandas.Series objects of equal length
    
    
        Returns
        ----------
        crosscorr : float
        """
        if wrap:
            shiftedy = datay.shift(lag)
            shiftedy.iloc[:lag] = datay.iloc[-lag:].values
            return datax.corr(shiftedy)
        else: 
            return datax.corr(datay.shift(lag))
    
    
    d1 = df['S1_Joy']
    d2 = df['S2_Joy']
    seconds = 5
    fps = 30
    rs = [crosscorr(d1,d2, lag) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
    offset = np.ceil(len(rs)/2)-np.argmax(rs)
    f,ax=plt.subplots(figsize=(14,3))
    ax.plot(rs)
    ax.axvline(np.ceil(len(rs)/2),color='k',linestyle='--',label='Center')
    ax.axvline(np.argmax(rs),color='r',linestyle='--',label='Peak synchrony')
    ax.set(title=f'Offset = {offset} frames\nS1 leads <> S2 leads',ylim=[.1,.31],xlim=[0,300], xlabel='Offset',ylabel='Pearson r')
    ax.set_xticklabels([int(item-150) for item in ax.get_xticks()])
    plt.legend()
    plt.show()
    

    上图中,我们可以从负坐标推断出,Subject 1(S1)信号在引导信号间的相互作用(当 S2 被推进了 47 帧的时候相关性最高)。但是,这个评估信号在全局层面会动态变化,例如在这三分钟内作为引导信号的信号就会如此。另一方面,我们认为信号之间的相互作用也许会波动得更加明显,信号是引导还是跟随,会随着时间而转换。

    为了评估粒度更细的动态变化,我们可以计算加窗的时间滞后互相关(WTLCC)。这个过程会在信号的多个时间窗内反复计算时间滞后互相关。然后我们可以分析每个窗口或者取窗口上的总和,来提供比较两者之间领导者跟随者互动性差异的评分。

    # 加窗的时间滞后互相关
    seconds = 5
    fps = 30
    no_splits = 20
    samples_per_split = df.shape[0]/no_splits
    rss=[]
    for t in range(0, no_splits):
        d1 = df['S1_Joy'].loc[(t)*samples_per_split:(t+1)*samples_per_split]
        d2 = df['S2_Joy'].loc[(t)*samples_per_split:(t+1)*samples_per_split]
        rs = [crosscorr(d1,d2, lag) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
        rss.append(rs)
    rss = pd.DataFrame(rss)
    f,ax = plt.subplots(figsize=(10,5))
    sns.heatmap(rss,cmap='RdBu_r',ax=ax)
    ax.set(title=f'Windowed Time Lagged Cross Correlation',xlim=[0,300], xlabel='Offset',ylabel='Window epochs')
    ax.set_xticklabels([int(item-150) for item in ax.get_xticks()]);
    
    
    # 滑动窗口时间滞后互相关
    seconds = 5
    fps = 30
    window_size = 300 #样本
    t_start = 0
    t_end = t_start + window_size
    step_size = 30
    rss=[]
    while t_end < 5400:
        d1 = df['S1_Joy'].iloc[t_start:t_end]
        d2 = df['S2_Joy'].iloc[t_start:t_end]
        rs = [crosscorr(d1,d2, lag, wrap=False) for lag in range(-int(seconds*fps-1),int(seconds*fps))]
        rss.append(rs)
        t_start = t_start + step_size
        t_end = t_end + step_size
    rss = pd.DataFrame(rss)
    
    
    f,ax = plt.subplots(figsize=(10,10))
    sns.heatmap(rss,cmap='RdBu_r',ax=ax)
    ax.set(title=f'Rolling Windowed Time Lagged Cross Correlation',xlim=[0,300], xlabel='Offset',ylabel='Epochs')
    ax.set_xticklabels([int(item-150) for item in ax.get_xticks()])
    plt.show()
    

    如上图所示,是将时间序列分割成了 20 个等长的时间段,然后计算每个时间窗口的互相关。这给了我们更细粒度的视角来观察信号的相互作用。例如,在第一个窗口内(第一行),右侧的红色峰值告诉我们 S2 开始的时候在引导相互作用。但是,在第三或者第四窗口(行),我们可以发现 S1 开始更多的引导相互作用。我们也可以继续计算下去,那么就可以得出下图这样平滑的图像。

    时间滞后互相关和加窗时间滞后互相关是查看两信号之间更细粒度动态相互作用的很好的方法,例如引导-追随关系以及它们如何随时间改变。但是,这样的对信号的计算的前提是假设事件是同时发生的,并且具有相似的长度,这些内容将会在下一部分涵盖。


    3. 动态时间扭曲 —— 同步长度不同的信号


    动态时间扭曲(DTW)是一种计算两信号间路径的方法,它能最小化两信号之间的距离。这种方法最大的优势就是他能处理不同长度的信号。最初它是为了进行语言分析而被发明出来,DTW 通过计算每一帧对于其他所有帧的欧几里得距离,计算出能匹配两个信号的最小距离。一个缺点就是它无法处理缺失值,所以如果你的数据点有缺失,你需要提前插入数据。

    为了计算 DTW,我们将会使用 Python 的 dtw 包,它将能够加速运算。

    from dtw import dtw,accelerated_dtw
    
    
    d1 = df['S1_Joy'].interpolate().values
    d2 = df['S2_Joy'].interpolate().values
    d, cost_matrix, acc_cost_matrix, path = accelerated_dtw(d1,d2, dist='euclidean')
    
    
    plt.imshow(acc_cost_matrix.T, origin='lower', cmap='gray', interpolation='nearest')
    plt.plot(path[0], path[1], 'w')
    plt.xlabel('Subject1')
    plt.ylabel('Subject2')
    plt.title(f'DTW Minimum Path with minimum distance: {np.round(d,2)}')
    plt.show()
    

    如图所示我们可以看到白色凸形线绘制出的最短距离。换句话说,较早的 Subject2 数据和较晚的 Subject1 数据的同步性能够匹配。最短路径代价是d=.33,可以用来和其他信号的该值做比较。


    4. 瞬时相位同步


    最后,如果你有一段时间序列数据,你认为它可能有振荡特性(例如 EEG 和 fMRI),此时你也可以测量瞬时相位同步。它也可以计算两个信号间每一时刻的同步性。这个结果可能会因人而异因为你需要过滤数据以获得你感兴趣的波长信号,但是你可能只有未经实践的某些原因来确定这些波段。为了计算相位同步性,我们需要提取信号的相位,这可以通过使用希尔伯特变换来完成,希尔波特变换会将信号的相位和能量拆分开(你可以在这里学习更多关于希尔伯特变换的知识)。这让我们能够评估两个信号是否同相位(两个信号一起增强或减弱)。

    from scipy.signal import hilbert, butter, filtfilt
    from scipy.fftpack import fft,fftfreq,rfft,irfft,ifft
    import numpy as np
    import seaborn as sns
    import pandas as pd
    import scipy.stats as stats
    def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        b, a = butter(order, [low, high], btype='band')
        return b, a
    
    
    
    
    def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        y = filtfilt(b, a, data)
        return y
    
    
    lowcut  = .01
    highcut = .5
    fs = 30.
    order = 1
    d1 = df['S1_Joy'].interpolate().values
    d2 = df['S2_Joy'].interpolate().values
    y1 = butter_bandpass_filter(d1,lowcut=lowcut,highcut=highcut,fs=fs,order=order)
    y2 = butter_bandpass_filter(d2,lowcut=lowcut,highcut=highcut,fs=fs,order=order)
    
    
    al1 = np.angle(hilbert(y1),deg=False)
    al2 = np.angle(hilbert(y2),deg=False)
    phase_synchrony = 1-np.sin(np.abs(al1-al2)/2)
    N = len(al1)
    
    
    # 绘制结果
    f,ax = plt.subplots(3,1,figsize=(14,7),sharex=True)
    ax[0].plot(y1,color='r',label='y1')
    ax[0].plot(y2,color='b',label='y2')
    ax[0].legend(bbox_to_anchor=(0., 1.02, 1., .102),ncol=2)
    ax[0].set(xlim=[0,N], title='Filtered Timeseries Data')
    ax[1].plot(al1,color='r')
    ax[1].plot(al2,color='b')
    ax[1].set(ylabel='Angle',title='Angle at each Timepoint',xlim=[0,N])
    phase_synchrony = 1-np.sin(np.abs(al1-al2)/2)
    ax[2].plot(phase_synchrony)
    ax[2].set(ylim=[0,1.1],xlim=[0,N],title='Instantaneous Phase Synchrony',xlabel='Time',ylabel='Phase Synchrony')
    plt.tight_layout()
    plt.show()
    

    瞬时相位同步测算是计算两个信号每一刻同步性的很好的方法,并且它不需要我们像计算滑动窗口相关性那样任意规定窗口宽度。

    展开全文
    FrankieHello 2021-09-13 00:35:25
  • 289KB u010197810 2019-08-09 08:43:35
  • 4星
    4KB lipyconstance 2012-08-21 03:21:04
  • 1.4MB weixin_38609128 2021-06-16 00:48:41
  • 为了检测时间序列的相关性,我们经常使用自相关互相关或归一化互相关互相关(Cross-Correlation) 互相关是两个不同时间序列的比较,以检测具有相同最大值和最小值的指标之间是否存在相关性。例如:“两个音频...

    为了检测时间序列的相关性,我们经常使用自相关,互相关或归一化互相关。

    互相关(Cross-Correlation)

    互相关是两个不同时间序列的比较,以检测具有相同最大值和最小值的指标之间是否存在相关性。例如:“两个音频信号同相吗?”

    为了检测两个信号之间的相关程度,我们使用互相关。 只需将两个时间序列相乘和相加即可计算得出。

    在以下示例中,序列A和B是互相关的,但序列C都不与此相关。

    cross-correlation

    a = [1, 2, -2, 4, 2, 3, 1, 0]
    b = [2, 3, -2, 3, 2, 4, 1, -1]
    c = [-2, 0, 4, 0, 1, 1, 0, -2]
    
    pyplot.figure()
    pyplot.title("Series")
    pyplot.plot(range(len(a)), a, color="#f44e2e", label="a")
    pyplot.plot(range(len(b)), b, color="#27ccc0", label="b")
    pyplot.plot(range(len(c)), c, color="#273ecc", label="c")
    pyplot.ylabel("value")
    pyplot.xlabel("index")
    pyplot.show()
    

    c r o s s _ c o r r ( x , y ) = ∑ i = 0 n − 1 x i ∗ y i cross\_corr(x, y) = \sum_{i=0}^{n-1} x_i*y_i cross_corr(x,y)=i=0n1xiyi

    使用上面的互相关公式,我们可以计算序列之间的相关度。

    import numpy as np
    def cross_corr(set1, set2):
        return np.sum(set1 * set2)
    
    a = np.array([1, 2, -2, 4, 2, 3, 1, 0])
    b = np.array([2, 3, -2, 3, 2, 4, 1, -1])
    c = np.array([-2, 0, 4, 0, 1, 1, 0, -2])
    
    print(f"Cross Correlation a,b: {cross_corr(a, b)}")
    print(f"Cross Correlation a,c: {cross_corr(a, c)}")
    print(f"Cross Correlation b,c: {cross_corr(b, c)}")
    
    OUTPUT:
    Cross Correlation a,b: 41
    Cross Correlation a,c: -5
    Cross Correlation b,c: -4
    

    序列a和b互相关性较大,值为41。序列a和c,a和b相关性较小,值分别为:-5,-4.

    标准化的互相关(Normalized Cross-Correlation)

    标准化的互相关也是用于两个时间序列的比较,但是使用了不同的评分结果。除了简单的互相关,它还可以比较具有不同值范围的指标。例如:“商店中的顾客数量与每天的销售数量之间是否存在相关性?”

    使用互相关值评价相关性存在三个问题:

    1. 互相关的评分值难以理解。
    2. 两个序列必须具有相同的幅度。如如果一个序列值缩小了一半,那么他的相关性会减少。 c r o s s _ c o r r ( a , a / 2 ) = 19.5 cross\_corr(a,a/2)=19.5 cross_corr(a,a/2)=19.5
    3. 无法解决序列值为0的问题,根据互相关公式,由于 0 ∗ 0 = 0 0*0=0 00=0 0 ∗ 200 = 0 0*200=0 0200=0,因此无法解决0值。

    为了解决这些问题,使用标准化的互相关

    n o r m _ c o r r ( x , y ) = ∑ n = 0 n − 1 x i ∗ y i ∑ n = 0 n − 1 x i 2 ∗ ∑ n = 0 n − 1 y i 2 norm\_corr(x,y)=\dfrac{\sum_{n=0}^{n-1} x_i*y_i}{\sqrt{\sum_{n=0}^{n-1} x_i^2 * \sum_{n=0}^{n-1} y_i^2}} norm_corr(x,y)=n=0n1xi2n=0n1yi2 n=0n1xiyi

    计算上边a,b,c序列的标准化的相关性如下:

    import numpy as np
    def norm_cross_corr(set1, set2):
        # python中求Normalized Cross Correlation的函数是: statsmodels.tsa.stattools.ccf
        return np.sum(set1 * set2) / (np.linalg.norm(set1) * np.linalg.norm(set2))
    
    a = np.array([1, 2, -2, 4, 2, 3, 1, 0])
    b = np.array([2, 3, -2, 3, 2, 4, 1, -1])
    c = np.array([-2, 0, 4, 0, 1, 1, 0, -2])
    
    print(f"Normalized Cross Correlation a,b: {norm_cross_corr(a, b)}")
    print(f"Normalized Cross Correlation a,c: {norm_cross_corr(a, c)}")
    print(f"Normalized Cross Correlation b,c: {norm_cross_corr(b, c)}")
    
    OUTPUT:
    Normalized Cross Correlation a,b: 0.9476128352180998
    Normalized Cross Correlation a,c: -0.15701857325533194
    Normalized Cross Correlation b,c: -0.11322770341445959
    

    标准化的互相关值很容易理解:

    1. 值越高,相关性越高。
    2. 当两个信号完全相同时,最大值为1: n o r m _ c r o s s _ c o r r ( a , a ) = 1 norm\_cross\_corr(a,a)=1 norm_cross_corr(a,a)=1
    3. 当两个信号完全相反时,最小值为-1: n o r m _ c r o s s _ c o r r ( a , − a ) = − 1 norm\_cross\_corr(a,-a)=-1 norm_cross_corr(a,a)=1

    标准化的互相关可以检测两个幅度不同的信号的相关性,如: n o r m _ c r o s s _ c o r r ( a , a / 2 ) = 1 norm\_cross\_corr(a,a/2)=1 norm_cross_corr(a,a/2)=1。信号A和具有一半振幅的同一信号之间具有最高的相关性!

    自相关(Auto-Correlation)

    自相关是时间序列在不同时间与其自身的比较。例如,其目的是检测重复模式或季节性。例如:“服务器网站上是否有每周的季节性?”“本周的数据与前一周的数据高度相关吗?”

    自相关的用途非常广泛,其中之一就是检测由于季节性导致的可重复的模式。

    下图展示了具有季节性为8的序列:

    image

    import numpy as np
    np.random.seed(10)
    a = np.tile(np.array(range(8)), 8) + np.random.normal(loc=0.0, scale=0.5, size=64)
    pyplot.figure(figsize=(12, 4))
    pyplot.title("Series")
    pyplot.plot(range(len(a)), a, color="#f44e2e", label="a")
    pyplot.ylabel("value")
    pyplot.xlabel("index")
    pyplot.legend(loc="upper right")
    pyplot.show()
    

    图像由如上代码生成,他是1~8的可重复序列,并且夹杂了一些随机噪声。

    接下来计算时间偏移为4和时间偏移为8时信号与自身之间的自相关。
    image

    import numpy as np
    np.random.seed(10)
    a = np.tile(np.array(range(8)), 8) + np.random.normal(loc=0.0, scale=0.5, size=64)
    
    ar4 = a[:len(a) - 4]
    ar4_shift = a[4:]
    print(f"Auto Correlation with shift 4: {cross_correlation(ar4, ar4_shift)}")
    
    pyplot.figure(figsize=(12, 4))
    pyplot.title("Series")
    pyplot.plot(range(len(ar4)), ar4, color="blue", label="ar4")
    pyplot.plot(range(len(ar4_shift)), ar4_shift, color="green", label="ar4_shift")
    pyplot.ylabel("value")
    pyplot.xlabel("index")
    pyplot.legend(loc="upper right")
    pyplot.show()
    
    OUTPUT:
    Auto Correlation with shift 4: 623.797612892277
    

    image

    import numpy as np
    np.random.seed(10)
    a = np.tile(np.array(range(8)), 8) + np.random.normal(loc=0.0, scale=0.5, size=64)
    
    ar8 = a[:len(a) - 8]
    ar8_shift = a[8:]
    print(f"Auto Correlation with shift 4: {cross_correlation(ar8, ar8_shift)}")
    
    pyplot.figure(figsize=(12, 4))
    pyplot.title("Series")
    pyplot.plot(range(len(ar4)), ar4, color="blue", label="ar4")
    pyplot.plot(range(len(ar4_shift)), ar4_shift, color="green", label="ar4_shift")
    pyplot.ylabel("value")
    pyplot.xlabel("index")
    pyplot.legend(loc="upper right")
    pyplot.show()
    
    OUTPUT:
    Auto Correlation with shift 8: 996.8417240253186
    

    上图清楚地显示了时间偏移8处的高自相关,而时间偏移4处没有。计算出的自相关值可以看到shift为8时的自相关值大于shift为4时的值。所以相比于4,我们觉得季节性更可能为8。

    标准化的自相关(Normalized Auto-Correlation)

    如前所述,用标准化的互相关可以更好地表达两个序列的相关性

    print(f"Auto Correlation with shift 4: {norm_cross_correlation(ar4, ar4_shift)}")
    print(f"Auto Correlation with shift 8: {norm_cross_correlation(ar8, ar8_shift)}")
    
    OUTPUT:
    Auto Correlation with shift 4: 0.5779124931215941
    Auto Correlation with shift 8: 0.9874907451139486
    

    相关性与时移(Correlation with Time Shift)

    标准化的互相关与时移(Normalized Cross-Correlation with Time Shift)

    时移可以应用于所有上述算法。想法是将指标与具有不同“时间偏移”的另一个指标进行比较。将时间偏移应用于归一化互相关函数将导致“具有X的时间偏移的归一化互相关”。这可以用来回答以下问题:“当许多顾客进我的商店时,我的销售额在20分钟后会增加吗?”

    为了处理时移,我们将原始信号与由x元素向右或向左移动的另一个信号关联。就像我们为自相关所做的一样。

    为了检测两个指标是否与时移相关,我们需要计算所有可能的时移。

    image

    image

    from statsmodels.tsa.stattools import ccf
    a = np.array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4])
    b = np.array([1, 2, 3, 3, 0, 1, 2, 3, 4, 0, 1, 1, 4, 4, 0, 1, 2, 3, 4, 0])
    res_ary = ccf(a, b, unbiased=False)  # 计算各种time shift下的互相关值
    
    print(res_ary)
    
    pyplot.bar(range(len(res_ary)), res_ary, fc="blue")
    pyplot.show()
    
    OUTPUT:
    [-3.29180823e-17  8.67261700e-01  1.27247799e-01 -4.65751654e-01
     -3.92862138e-01 -1.09726941e-17  6.20178595e-01  1.27247799e-01
     -2.92793480e-01 -2.94028896e-01 -2.47083106e-02  3.48387179e-01
      1.02539489e-01 -1.19835306e-01 -1.70487343e-01 -2.47083106e-02
      1.01304073e-01  5.31228677e-02 -2.10020640e-02 -4.69457901e-02]
    

    从上图以及数据可以看出,当time shift为1时,序列a和b高度相关。

    标准化的自相关与时移(Normalized Auto-Correlation with Time Shift)

    image

    image

    from statsmodels.tsa.stattools import acf
    
    np.random.seed(5)
    a = np.tile(np.array(range(8)), 8) + np.random.normal(loc=0.0, scale=0.5, size=64)
    res_ary = acf(a, nlags=30, fft=False)
    print(res_ary)
    pyplot.bar(range(len(res_ary)), res_ary, fc="blue")
    pyplot.show()
    
    OUTPUT:
    [ 1.          0.33436187 -0.06325233 -0.39112367 -0.46168766 -0.43537856
     -0.15175395  0.24403119  0.86075234  0.31051923 -0.03270547 -0.32083613
     -0.40031665 -0.38383002 -0.14645852  0.20649692  0.72864895  0.2693097
     -0.01707009 -0.2658133  -0.3339456  -0.32771517 -0.1214095   0.17008013
      0.60740024  0.22492556 -0.01633838 -0.21157885 -0.27896827 -0.27765971
     -0.10762048]
    

    原始序列是有8的季节性,自相关检测也看到了在8的整数倍的偏移处得到了较高的自相关值。

    将相关指标聚类(Cluster Correlated Metrics Together)

    我们可以用Normalized Cross Correlation,根据相似度将相关性较大的metric聚类。

    使用数据集:graphs45.csv

    import numpy as np
    import pandas as pd
    
    def get_correlation_table(metric_df):
        metric_cnt = metric_df.shape[1]
        correlation_table = np.zeros((metric_cnt, metric_cnt))
        for i in range(metric_cnt):
            metric_1 = metric_df.iloc[:, i]
            for j in range(metric_cnt):
                if i == j:
                    continue
                else:
                    metric_2 = metric_df.iloc[:, j]
                    cc_ary = ccf(metric_1, metric_2, unbiased=False)
                    correlation_table[i, j] = cc_ary[0]
        return correlation_table
    
    
    def find_related_metric(correlation_table, orig, high_corr):
        metric_corr_table = correlation_table[orig]
        corr_metric_lst = []
        for i in range(len(metric_corr_table)):
            if metric_corr_table[i] > high_corr:
                corr_metric_lst.append(i)
    
        return corr_metric_lst
    
    if __mian__:
        metric_df = pd.read_csv("./graphs45.csv")
        correlation_table_ary = get_correlation_table(metric_df)
    
        orig = 3
        high_corr = 0.9
        corr_metric_lst = find_related_metric(correlation_table_ary, orig, high_corr)
        corr_metric_lst.append(orig)
        print(corr_metric_lst)
    
        pyplot.figure()
        pyplot.title("Series")
        for idx in corr_metric_lst:
            metric = metric_df.iloc[:, idx]
            pyplot.plot(range(len(metric)), metric, label=f"graph_{idx + 1}")
        pyplot.ylabel("value")
        pyplot.xlabel("index")
        pyplot.legend(loc="upper right")
        pyplot.show()
    

    get_correlation_table()计算了所有metric之间的标准化的互相关值,结果存在correlation_table中,correlation_table[i,j]代表第i个metric和第j个metric之间的相关值。

    然后find_related_metric()根据相关值表找出和graph_4相关值大于0.9的序列,最后找出的相关序列画图如下:

    image

    参考

    1. Understanding Cross-Correlation, Auto-Correlation, Normalization and Time Shift
    2. Detecting Correlation Among Multiple Time Series
    展开全文
    Eaton18 2019-11-18 00:52:05
  • r软件时间序列分析论文 数据科学 , 机器学习 (Data Science, Machine Learning) In machine learning with time series, using features extracted from series is more powerful than simply treating a time ...

    r软件时间序列分析论文

    数据科学机器学习 (Data Science, Machine Learning)

    In machine learning with time series, using features extracted from series is more powerful than simply treating a time series in a tabular form, with each date/timestamp in a separate column. Such features can capture the characteristics of series, such as trend and autocorrelations.

    在具有时间序列的机器学习中,使用从序列中提取的特征比仅以表格形式处理时间序列(每个日期/时间戳在单独的列中)更强大。 这些特征可以捕获序列的特征,例如趋势和自相关。

    But… what sorts of features can you extract and how do you select among them?

    但是……您可以提取哪些类型的特征,以及如何在其中进行选择?

    In this article, I discuss the findings of two papers that analyze feature-based representations of time series. The papers conduct comprehensive work to collect thousands of time series feature extractors and evaluate which features capture the most useful information from a series.

    在本文中,我讨论了两篇分析基于特征的时间序列表示的论文的发现。 这些论文进行了全面的工作,以收集成千上万个时间序列特征提取器,并评估哪些特征捕获了序列中最有用的信息。

    The papers show how to compare time series by extracting features that describe the series behavior and suggest a pipeline for identifying an “optimal” subset of time series features.

    这些论文展示了如何通过提取描述序列行为的特征并建议用于识别时间序列特征的“最佳”子集的管道来比较时间序列。

    为什么这很重要? (Why Is This Important?)

    There are two basic ways to compare time series:

    有两种比较时间序列的基本方法:

    1. A similarity measure that quantifies whether two-time series are close (on average) across time, such as Dynamic Time Warping. These measures are typically best for short, aligned series of equal length. They tend to have poor scalability, with quadratic computation in both the number of time series and series length because distances must be computed between all pairs.

      一种用于量化两个时间序列在整个时间上是否接近(平均)的相似性度量 ,例如Dynamic Time Warping 。 这些措施通常最适合短而对齐的等长序列。 它们往往具有较差的可伸缩性,因为在时间序列的数量和序列长度上都需要进行二次计算,因为必须在所有对之间计算距离。

    2. Define similarity between series in terms of features extracted from time series using time series analysis algorithms. Feature extractors do not require series to be of equal length. The result is an interpretable summary of the dynamical characteristics of each series. These features can then be used for machine learning.

      使用时间序列分析算法从时间序列提取特征方面定义序列之间的相似性。 特征提取器不需要序列的长度相等。 结果是每个系列动力学特性的可解释性总结。 这些功能可以用于机器学习。

    Interpretability is another key: time series features can capture complex, time-varying patterns in a set of interpretable characteristics.

    可解释性是另一个关键:时间序列特征可以以一组可解释的特征捕获复杂的时变模式。

    Problematically, there are a vast number of methods to extract interpretable features from time series. Further, feature selection is often done manually and subjectively.

    有问题的是,有很多方法可以从时间序列中提取可解释的特征。 此外,特征选择通常是手动和主观地完成的。

    What sort of features can be extracted from series and how could you select among them?

    可以从系列中提取什么样的特征,如何从中选择?

    Image for post
    hctsa: A Computational Framework for Automated Time-Series Phenotyping Using Massive Feature Extraction hctsa:使用大规模特征提取进行自动时间序列表型分析的计算框架

    高度比较的时间序列分析:时间序列的经验结构及其方法 (Highly comparative time-series analysis: the empirical structure of time series and their methods)

    Paper motivation: although time series are studied across scientific disciplines (e.g. stock prices in finance, human heartbeats in medicine), different methods for time series analysis have been developed separately in different disciplines.

    论文动机:尽管跨学科研究了时间序列(例如金融中的股票价格,医学上的人的心跳),但在不同学科中分别开发了不同的时间序列分析方法

    Given the great number of methods, it is difficult to determine how methods developed by different disciplines are related. As a result, how can a practitioner select the optimal method for their data?

    鉴于方法众多,因此很难确定不同学科开发的方法之间的关系。 结果,从业者如何为他们的数据选择最佳方法?

    To address this challenge, the HCTSA paper…

    为了应对这一挑战,HCTSA论文…

    • Assembles an extensive annotated library of time series data and methods for time series analysis.

      组装了一个广泛的带注释的时间序列数据库和时间序列分析方法。
    • Models time series methods according to their behavior on the data and group time series by their measured properties.

      根据时间序列方法在数据上的行为对时间序列方法进行建模,并通过其测量属性将时间序列分组。
    • Introduces a range of comparative analysis techniques for series and their methods. First, the ability to link given time series to similar real-world and model-generated series. Second, the ability to link specific time series analysis methods to a range of alternatives across the literature.

      介绍了一系列用于系列及其方法的比较分析技术。 首先,可以将给定的时间序列链接到类似的真实世界和模型生成的序列。 其次,将特定的时间序列分析方法链接到整个文献中的其他方法的能力。

    HCTSA框架和范围 (HCTSA Framework and Scope)

    The paper is scope extensive: the authors annotated a library of 38,190 univariate time series and 9,613 time series analysis algorithms.

    本文涉及面很广:作者注释了38,190个单变量时间序列和9,613个时间序列分析算法的库。

    The time series analysis methods vary in form, ranging from summary statistics to statistical model fits. Each transformation summarizes an input series with a single real number.

    时间序列分析方法的形式各不相同,从汇总统计信息到统计模型拟合不等。 每个转换都会汇总一个具有单个实数的输入序列。

    The library of time series transformations cover a wide range of time series properties:

    时间序列转换库涵盖了广泛的时间序列属性:

    • basic statistics of the distribution (e.g. location, spread, outlier properties)

      分布的基本统计信息(例如位置,分布,离群值属性)
    • linear correlations (e.g. autocorrelations, features of the power spectrum)

      线性相关(例如,自相关,功率谱的特征)
    • stationarity (e.g. sliding window measures, unit root tests)

      平稳性(例如,滑动窗口度量,单位根检验)
    • information-theoretic and entropy measures (e.g. auto-mutual information, Approximate Entropy)

      信息理论和熵测度(例如,自动互信息,近似熵)
    • methods from the physical nonlinear time-series analysis literature (e.g. correlation dimension)

      物理非线性时间序列分析文献中的方法(例如,相关维)
    • linear and nonlinear model fits (e.g. goodness of fit and parameters from autoregressive models)

      线性和非线性模型拟合(例如拟合优度和自回归模型的参数)
    • others (e.g. wavelet methods)

      其他(例如小波方法)

    For transformations that require parameter values, the transformation is repeated for multiple parameters. A single “operation” is considered a transformation plus a single parameter value. Of the 9k operations evaluated in the paper, a single transformation might be counted multiple times, once for each parameter value. The paper evaluates approximately 1k unique transformations.

    对于需要参数值的转换,将对多个参数重复该转换。 单个“操作”被视为转换加上单个参数值。 在本文评估的9k运算中,单个转换可能会被计数多次,每个参数值一次。 本文评估了大约1k个唯一转换。

    HCTSA:时间序列分析方法的经验结构 (HCTSA: Empirical structure of time series analysis methods)

    Image for post
    Figure 3 A. A summary of the four main classes of time series operations.
    图3 A.时间序列操作的四个主要类别的摘要。

    The authors used k-medoids clustering to identify four broad categories of time series analysis operations:

    作者使用k-medoids聚类来识别时间序列分析操作的四大类:

    1. Linear correlation

      线性相关
    2. Stationarity (Properties that change with time)

      平稳性(随时间变化的属性)
    3. Information theory

      信息论
    4. Nonlinear time series analysis

      非线性时间序列分析

    The clustering analysis revealed that a subset of 200 time series operations, or an empirical fingerprint of a series’ behavior, can approximate the 8,651 operations considered. The 200 operations summarize different behaviors of time series analysis methods. These operations include techniques developed in a variety of disciplines.

    聚类分析表明,200个时间序列操作的子集或序列行为经验指纹可以近似考虑所考虑的8,651个操作。 这200个操作总结了时间序列分析方法的不同行为。 这些操作包括在各种学科中开发的技术。

    Further, the analysis uncovered a local structure surrounding each target operation. For a given operation, they were able to identify alternative operations with similar behavior.

    此外,分析发现了围绕每个目标操作的局部结构。 对于给定的操作,他们能够识别行为相似的替代操作。

    Image for post
    Figure 3 B. “A network representation of the operations in our library that are most similar to the Approximate Entropy algorithm, “ApEn( 2,0.2)”…, which were retrieved from our library automatically. Each node in the network represents an operation and links encode distances between them…. Annotated scatter plots show the outputs of Approximate Entropy against a representative member of each shaded community (indicated by a heavily outlined node).”
    图3 B.“我们库中与最近似熵算法“ ApEn(2,0.2)”最相似的操作的网络表示,…是从我们库中自动检索的。 网络中的每个节点代表一个操作,并且链接对它们之间的距离进行编码……。 带注释的散点图显示了每个阴影社区的代表成员的近似熵输出(由轮廓突出的节点指示)。”

    “By comparing their empirical behaviour, the techniques demonstrated above can be used to connect new methods to alternatives developed in other fields in a way that encourages interdisciplinary collaboration on the development of novel methods for time-series analysis that do not simply reproduce the behaviour of existing methods” [1]

    “通过比较他们的经验行为,上面展示的技术可以用于将新方法与其他领域开发的替代方法联系起来,从而鼓励跨学科合作,开发时间序列分析的新方法,而不仅仅是再现行为的行为。现有方法” [1]

    HCTSA:时间序列的经验结构 (HCTSA: Empirical structure of time series)

    Time series can be represented by properties that capture important dynamical behavior of the series. The authors use 200 representative operations to compare 24,577 time series from different systems and of varying lengths.

    时间序列可以由捕获序列的重要动力学行为的属性表示。 作者使用200个代表性操作来比较来自不同系统和不同长度的24,577个时间序列。

    This empirical fingerprint of 200 diverse time-series analysis operations facilitates a meaningful comparison of scientific time series.

    200种不同时间序列分析操作的经验指纹有助于对科学时间序列进行有意义的比较。

    To group their library of 24k time series, the authors used complete linkage clustering to form 2,000 clusters. Due to the wide range of time series properties used, the clusters grouped series according to dynamics, even when the lengths differ.

    为了将他们的24k时间序列库分组,作者使用了完整的链接聚类来形成2,000个聚类。 由于使用了广泛的时间序列属性,因此即使长度不同,聚类也会根据动力学将序列分组。

    Most clusters grouped time series measured from the same system:

    大多数群集将从同一系统测得的时间序列分组:

    Image for post
    Figure 4A: “Most clusters formed in this way are homogenous groups of time series of a given real-world or model system”
    图4A:“以这种方式形成的大多数集群是给定的现实世界或模型系统的时间序列的同质组”

    Some clusters contained series generated by different systems:

    一些集群包含由不同系统生成的序列:

    Image for post
    Figure 4B: “A time-series cluster is plotted that contains time series generated by three different iterative maps with parameters that specify a common recurrence relationship. Time-series segments of 150 samples are plotted and labeled with the parameter A of the map that generated them.”
    图4B:绘制了一个时间序列簇,其中包含由三个不同的迭代图生成的时间序列,这些迭代图的参数指定了共同的递归关系。 将绘制150个样本的时间序列片段,并用生成它们的图的参数A进行标记。”

    The reduced representation of time series allows you to retrieve a local neighborhood of series with similar properties. This allows you to automatically relate real-world time series to similar, model-generated time series.

    时间序列的简化表示使您可以检索具有相似属性的序列的局部邻域。 这使您可以自动将现实世界的时间序列与模型生成的类似时间序列相关联。

    Thus, the transformations can be used to suggest suitable families of models for use in real-world systems.

    因此,这些转换可用于建议适用于实际系统的模型族。

    Image for post
    Fig 4C. Opening share price series for Oxford Instruments (OXIG) in big red point; the most similar real world time series are opening share prices of other stocks (red nodes). Most similar model-generated TS are from stochastic differential equations (blue nodes. Links in network rep similarities between the series according to Euclidean distances between the normalized feature vectors.
    图4C。 牛津仪器(OXIG)的开盘股价序列具有较大的红点; 现实世界中最相似的时间序列是其他股票(红色节点)的开盘价。 大多数类似的模型生成的TS均来自随机微分方程(蓝色节点。根据归一化特征向量之间的欧几里得距离,序列之间的网络重复性相似性链接。

    HCTSA守则 (HCTSA Code)

    The code for Highly Comparative Time Series Analysis can be found on GitHub; however, it is written in Matlab. (You can use the hctsa package from python using the pyopy package). The hctsa package allows thousands of features to be extracted from a time series. The software also has an accompanying paper.

    可在GitHub上找到高度比较时间序列分析的代码; 但是,它是用Matlab编写的。 (您可以使用pyopy包从python使用pyopy包)。 hctsa包允许从一个时间序列中提取成千上万个功能。 该软件还附有论文

    Of important note, it is slow to run. Reducing the full set of HCTSA operations to even 200 of the thousands of candidate features is computationally expensive. This approach is infeasible for some applications, especially those with large training data.

    重要的是,它运行缓慢。 将全套HCTSA操作减少到数千个候选特征中的200个在计算上是昂贵的。 对于某些应用程序,尤其是具有大量训练数据的应用程序,这种方法是不可行的。

    HCTSA also has a web platform, CompEngine. CompEngine “is a self-organizing database of time-series data that allows users to upload, explore, and compare thousands of diverse types of time-series data.” [4]

    HCTSA还具有一个Web平台CompEngine 。 CompEngine“是一个时间序列数据的自组织数据库,允许用户上载,浏览和比较数千种不同类型的时间序列数据。” [4]

    Image for post
    Ewan Munro on Ewan Munroflickr flickr上的照片

    catch22,CAnonical时间序列特征 (catch22, CAnonical Time-series CHaracteristics)

    The subsequent catch22: CAnonical Time-series CHaracteristics paper (2019) builds on HCTSA by reducing the set of representative features to 22 time series features that:

    随后的内容22:CAnonical时间序列Characteristics论文(2019)建立在HCTSA的基础上, 将代表性特征的集合减少到22个时间序列特征 ,这些特征包括:

    1. exhibit strong classification performance across a given collection of time-series problems, and

      在给定的时间序列问题集合中表现出强大的分类性能,并且
    2. are minimally redundant, and

      最少冗余,并且
    3. capture the diversity of analysis contained in HCTSA.

      捕获HCTSA中包含的分析多样性。

    The paper creates a data-driven subset of the most useful features extracted from a time series. The authors compare across a diverse set of time series analysis algorithms, starting with the features in the HCTSA toolbox.

    本文创建了从时间序列中提取的最有用功能的数据驱动子集。 作者从HCTSA工具箱中的功能开始,对各种时间序列分析算法进行了比较。

    The catch22 time series characteristics capture a diverse and interpretable time series “signature” based on their properties.

    catch22时间序列特征基于其特性捕获了多种且可解释的时间序列“签名”。

    This signature includes linear and non-linear temporal auto-correlation, successive differences, value distributions and outliers, and fluctuation scaling properties.

    该签名包括线性和非线性时间自相关,连续差异,值分布和离群值以及波动比例属性。

    catch22功能的好处 (Benefits of catch22 features)

    • Fast computation (~1000x faster than full HCTSA feature set in Matlab)

      快速计算(比Matlab中完整的HCTSA功能集快1000倍)
    • Provides low dimensional summary of time series

      提供时间序列的低维摘要
    • Interpretable characteristics that are useful for classification and clustering.

      可解释的特征,对分类和聚类很有用。

    Further, if the catch22 features are not appropriate for your problem, the feature selection pipeline is general. The pipeline can be used to select informative subsets of features new or more complex problems.

    此外,如果catch22功能不适合您的问题,则功能选择管道很通用。 管道可用于选择新的或更复杂问题的特征性信息子集。

    Catch22功能评分 (Catch22 feature scoring)

    The authors score features by evaluating decision tree classification accuracy across a set of 93 classification problems from the Time Series Classification Repository. Performance with 4791 features from HCTSA has 77.2% mean class-balanced accuracy across all tasks. Performance with smaller set of 22 features is 71.7% mean class-balanced accuracy.

    作者通过评估时间序列分类库中的93个分类问题的决策树分类准确性来为特征评分。 HCTSA具有4791功能的性能在所有任务中具有77.2%的平均班级平衡准确性。 具有22个功能的较小集合的性能为71.7%的平均类平衡准确性。

    Catch22功能选择管道 (Catch22 feature selection pipeline)

    For all data sets, each time series feature was linearly rescaled to unit 0–1 interval. This scaling may not be appropriate for some real-world applications.

    对于所有数据集,每个时间序列特征均线性调整为单位0–1间隔。 这种缩放可能不适用于某些实际应用。

    First, the authors excluded features sensitive to mean and variance of distribution of values because the majority of series were normalized.

    首先,作者排除了对值的均值和方差敏感的特征,因为大多数序列都已归一化。

    For some applications, this preselection is not desirable. If working with non-normalized series, you should consider including the distributional features, such as mean and standard deviation. These can lead to significant performance gains.

    对于某些应用,这种预选择是不希望的。 如果使用非归一化序列,则应考虑包括分布特征,例如均值和标准差。 这些可以导致显着的性能提升。

    Next, the authors excluded the transformations that frequently output special values. Special values indicate that an algorithm is not suitable for the input data, or that it did not evaluate successfully.

    接下来,作者排除了经常输出特殊值的转换。 特殊值表示算法不适合输入数据,或者评估失败。

    Last, the authors created a pipeline to filter for features that can individually discriminate across a range of real-world data. The pipeline then filtered for those that have complementary behavior.

    最后,作者创建了一个管道,以筛选可分别区分一系列实际数据的功能。 然后,管道会筛选出具有互补行为的管道。

    The feature selection pipeline had 3 rounds:

    功能选择管道进行了三轮:

    1. Statistical pre-filtering: filter out features whose performance were statistically insignificant on the given learning tasks.

      统计预过滤:过滤掉在​​给定学习任务中性能在统计上不重要的特征。
    2. Performance filtering: select features that perform best across all tests. “Performance” is the ability to distinguish between labeled classes in 93 classification tasks with a decision tree classifier.

      性能过滤:选择在所有测试中性能最好的功能。 “性能”是使用决策树分类器区分93个分类任务中标记的类的能力。
    3. Redundancy minimization. The top features were clustered (hierarchical clustering with complete linkage) into groups according to performance scores across tasks. From each cluster, a single representative feature was selected for the feature set. The representative feature selected as the one with highest score across tasks — unless it was computationally intensive, in which case another high-accuracy feature with greater interpretability and efficiency was manually selected.

      冗余最小化。 根据任务之间的性能得分,将主要功能(通过完全链接的层次化群集)进行分组。 从每个群集中,为功能集选择一个代表性功能。 代表性特征被选为在所有任务中得分最高的特征-除非计算量大,否则将手动选择另一种具有更高可解释性和效率的高精度特征。
    Image for post
    Springer) 施普林格 )

    准确性/可解释性的权衡 (Accuracy / Interpretability Trade-off)

    The authors compared classification performance using the catch22 features with a wide variety of time series classification algorithms, such as those implemented in sktime.

    作者将使用catch22功能的分类性能与各种时间序列分类算法(例如在sktime实现的算法)进行了sktime

    The classification of time series with catch22 features, despite large dimensionality reduction, results in “similar” performance to alternative methods. The authors admit that majority of datasets exhibit better performance using existing algorithms than catch22.

    尽管具有较大的降维效果,但具有catch22特征的时间序列分类却导致与替代方法“相似”的性能。 作者承认, 使用现有算法大多数数据集表现出比catch22更好的性能。

    The paper often claims that catch22 only has a “small’ reduction in accuracy. (The authors did not publish the performance of classifiers with catch22 features). In one instance, they called a decrease from 99.2% to 89.5% “small”, but in my opinion, this is not small for many applications.

    该论文经常声称catch22的准确性仅“小”降低。 (作者未发布具有catch22功能的分类器的性能)。 在一种情况下,他们称从99.2%降低到89.5%是“小”,但在我看来,这对于许多应用程序来说并不小。

    While the authors failed to prove, in my view, that a classification model built with the catch22 features could outperform a native time series classifier, catch22 does offer interpretable features for model explanation.

    在我看来,尽管作者未能证明使用catch22特征构建的分类模型可以胜过本机时间序列分类器, 但是catch22确实提供了可解释的特征用于模型解释

    In particular, the authors highlighted one classifier where a single feature was able to perfectly separate two classes (series = triangle or noise). The feature “quantifies the length of the longest continued descending increments in the data”. Clearly, this is simple to explain.

    尤其是,作者强调了一个分类器,其中一个功能可以完美地将两个分类(序列=三角形或噪声)分开。 该功能“量化数据中最长的连续下降增量的长度”。 显然,这很容易解释。

    Image for post
    Fig 8A. Two classes (triangle vs noise series) are perfectly separable according to a single feature.
    图8A。 根据单个功能,两个类别(三角形与噪声系列)是完全可分离的。

    备用时间序列功能集 (Alternative Time Series Feature Sets)

    The authors noted that “There is no single representation that is best for all time-series datasets.” Instead, “the optimal representation depends on the structure of the dataset and the questions being asked of it.” [3]

    作者们指出:“没有一种最适合所有时间序列数据集的表示形式。” 相反,“最佳表示形式取决于数据集的结构和所要提出的问题。” [3]

    Thus, the catch22 features may not be the optimal features for all time series datasets and tasks.

    因此,catch22特征可能不是所有时间序列数据集和任务的最佳特征。

    The catch22 feature representation often outperforms datasets that do not have “reliable shape differences between classes” relative to classifiers based on time-domain distance metrics.

    catch22特征表示相对于基于时域距离度量的分类器,其性能通常优于没有“可靠的类间形状差异”的数据集。

    The authors compared performance of catch22 features to the time series features available in the tsfeatures R package. On the same set of classification tasks, tsfeatures features had a 69.4% mean accuracy, compared to catch22’s 71.7% accuracy.

    作者将catch22功能的性能与tsfeatures R软件包中可用的时间序列功能进行了tsfeatures 。 在同一组分类任务中, tsfeatures特征的平均准确度为69.4%,而catch22的平均准确度为71.7%。

    实作 (Implementation)

    Extraction of the catch22 features has been implemented in C, with wrappers in Python, R, Matlab. An open-source implementation of catch22 can be found on GitHub.

    catch22功能的提取已在C中实现,并在Python,R,Matlab中使用了包装器。 catch22的开源实现可以在GitHub上找到

    The C version of catch22 exhibits near-linear computational complexity, O(N1.16) for time series length. For a time series with 10,000 observations, the catch22 can be computed in 0.5 seconds.

    catch22的C版本显示时间序列长度的近似线性计算复杂度O(N1.16)。 对于具有10,000个观测值的时间序列,可以在0.5秒内计算catch22。

    The code for the feature selection pipeline that produced the 22 features is available on GitHub at https://github.com/chlubba/op_importance.

    GitHub上的https://github.com/chlubba/op_importance上提供了用于生成22个功能的功能选择管道的代码。

    Image for post
    Image by Ann H on Pexels
    Ann HPexels上的 图片

    适用于实际问题 (Application to real problems)

    A wide range of features can be extracted from time series that describe the many properties and dynamics of a series.

    可以从时间序列中提取各种各样的特征,这些特征描述了序列的许多特性和动力学。

    The features analyzed in the HCTSA paper and are available on GitHub are comprehensive and informative. The key challenge is that there are “too many” features for most applications.

    HCTSA论文中分析的功能以及可以在GitHub上获得的功能都是全面且信息丰富的。 关键的挑战是大多数应用程序的功能太多。

    The catch22 features are tailored to capture key properties of the UCR/UEA datasets, which are short and phase aligned. The feature selection method could be rerun to generate reduced feature sets tailored to other applications.

    catch22的功能经过定制,可以捕获UCR / UEA数据集的关键属性,这些属性很短且相位对齐。 可以重新运行功能选择方法以生成适合其他应用程序的精简功能集。

    Indeed, new feature selection may be necessary in many applications where the series have different properties, such as those where location and variance of a data distribution are highly relevant. Distributional features were excluded from the catch22 analysis because the data considered were normalized. (Normalization removes location and shift).

    确实,在一系列具有不同属性的应用程序中,例如在数据分布的位置和方差高度相关的那些应用程序中,可能需要新的特征选择。 catch22分析排除了分布特征,因为考虑的数据已标准化。 (归一化删除位置和移位)。

    最后的话 (A Final Word)

    If you enjoyed this article, please follow me for more content about time series machine learning. Articles on time series classification and a taxonomy of time series features are in the works.

    如果您喜欢本文,请关注我以获取有关时间序列机器学习的更多内容。 有关时间序列分类和时间序列特征分类的文章正在撰写中。

    翻译自: https://medium.com/towards-artificial-intelligence/highly-comparative-time-series-analysis-a-paper-review-5b51d14a291c

    r软件时间序列分析论文

    展开全文
    weixin_26746401 2020-09-27 02:28:24
  • qtlyx 2016-12-05 16:27:12
  • vinking9393 2016-06-07 18:26:44
  • m0_46510245 2021-08-05 09:12:44
  • weixin_42294495 2021-04-21 09:46:30
  • m0_49271156 2021-03-15 17:47:30
  • m0_37876745 2021-04-30 11:20:36
  • tonyhxhd 2017-11-16 03:04:50
  • weixin_30588655 2021-04-20 04:47:02
  • fangsfdavid 2021-07-16 15:27:34
  • 952KB weixin_38750861 2021-01-27 17:39:13
  • 178KB weixin_38744270 2019-09-20 22:23:42
  • weixin_39645308 2021-05-23 11:07:38
  • a321123b 2020-12-19 21:43:55
  • weixin_42711949 2020-08-04 18:05:17
  • weixin_45577825 2021-09-17 21:37:42

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 19,609
精华内容 7,843
关键字:

时间序列互相关分析