Fibonacci 序列

在开始讲解时间序列的自回归模型（AutoRegression Model）之前，我们需要先介绍一下线性代数的基础知识。为了介绍线性代数的基础知识，我们可以先看一个简单的例子。考虑整数序列 Fibonacci 序列，它的通项公式是一个递归函数，每项的取值是由前两项所生成的，其具体的公式就是

，

其初始值是 。按照其递归公式来计算，我们可以详细写出前面的几项，那就是：

0，1，1，2，3，5，8，13，21，34，55，89，144，…

但是，计算 Fibonacci 的通项公式则比计算等差数列和等比数列的通项公式复杂的多，因为这里需要使用线性代数的技巧才能解决这个问题。

求解 Fibonacci 序列的通项公式 －－－ 矩阵对角化

根据 Fibonacci 数列的递归公式，基于矩阵乘法的定义，Fibonacci 序列可以写成如下形式：

这就意味着我们需要计算出矩阵 A 的幂。在线性代数里面，为了计算矩阵的 n 次方，除了通过矩阵乘法的公式直接计算之外，还有一个经典的技巧，那就是将矩阵对角化。详细来说，如果  的矩阵 A 能够对角化，那就是存在可逆矩阵 P 使得

其中  表示一个  的对角矩阵，其对角的元素从左上角到右下角依次是 。如果把矩阵 P 写成列向量的形式，i.e. ，那么以上的矩阵方程就可以转换为 , 。进一步来说，如果要计算矩阵 A 的幂，就可以得到：

另外，如果想知道一个矩阵 A 的特征值和特征向量，则需要计算以下多项式的根，i.e. 计算关于  的多项式的解，

其中 I 是一个单位矩阵（identity matrix）。

按照以上的思路，如果令

,

可以计算出 A 的两个特征值分别是  和 ，它们所对应的特征向量分别是：

.

因此直接计算可以得到

.

通过上面的计算方法，为了计算 Fibonacci 数列的通项公式，我们可以先把它转换成一个矩阵求幂的问题，于是我们就能矩阵对角化的方法把 Fibonacci 数列的通项公式求出来。

时间序列的弱平稳性

要讲解自回归模型，就必须提到时间序列的弱平稳性。一个时间序列  具有弱平稳性（Weak Stationary）指的是：

1.  对于所有的  都是恒定的；
2.  对于所有的  都是恒定的；
3.  与  的协方差对于所有的  都是恒定的。

另外，时间序列的自相关方程（AutoCorrelation Function）指的是对于 ，可以定义 ACF 为

.

如果时间序列  在弱平稳性的假定下，ACF 将会简化为

.

时间序列的自回归模型（AutoRegression Model）

AR(1) 模型

AR(1) 模型指的是时间序列  在时间戳  时刻的取值  与时间戳  时刻的取值  相关，其公式就是：

，

这个时间序列  满足如下条件：

1. ，并且  满足 iid 条件。其中  表示 Gauss 正态分布，它的均值是0，方差是 。
2.  与  是相互独立的（independent）。
3.  是弱平稳的，i.e. 必须满足 。

如果选择初始条件 ，则可以得到一些 AR(1) 模型的例子如下图所示：

从 AR(1) 以上的定义出发，我们可以得到：

1. .
2. .
3. .

Proof of 1. 从 AR(1) 的模型出发，可以得到

,

从而，.

Proof of 2. 从 AR(1) 的模型出发，可以得到

,

从而，.

Proof of 3. 令 . 从  的定义出发，可以得到：

从而，

.

AR(1) 模型与一维动力系统

特别的，如果假设  恒等于零，就可以得到  对于所有的  都成立。也就是可以写成一个一维函数的迭代公式：

下面我们要计算  的收敛性，这里的  表示函数  的  次迭代。

Method 1. 

通过  的定义直接计算可以得到：

,

令 ，可以得到 。这与  其实是保持一致的。

另外，如果 ，可以从公式上得到  当 。

Method 2.

将原函数转换成 Lipschitz 函数的形式，i.e. 如果 ，那么

.

因为 ，我们可以得到 . i.e.  趋近于 .

反之，如果 ，很容易得到

.

因此，在  这种条件下， as . 特别地，对于一阶差分方程  而言，如果 ，那么  的取值会越来越大，这与现实的状况不相符，所以在时间序列的研究中，一般都会假设 。

AR(p) 模型

按照之前类似的定义，可以把 AR(1) 模型扩充到 AR(p) 模型，也就是说：

1. AR(1) 模型形如：

2. AR(2) 模型形如：

3. AR(p) 模型形如：

AR(p) 模型的稳定性 －－－ 基于线性代数

对于 AR(2) 模型，可以假定  并且忽略误差项，因此可以得到简化版的模型形如：

.

写成矩阵的形式形如：

求解其特征多项式则是基于 ，求解可以得到 ，i.e. 。当  都在单位圆内部的时候，也就是该模型  满足稳定性的条件。

对于更加一般的 AR(p) 模型，也就是考虑 p 阶差分方程

可以用同样的方法将其转换成矩阵的形式，那就是：

计算 ，可以得到其特征多项式为：

当每个特征值都在单位圆盘内部的时候，i.e. ，该 p 阶差分方程

存在稳定性的解。

在之前的时间序列相似度算法中，时间戳都是一一对应的，但是在实际的场景中，时间戳有可能出现一定的偏移，但是两条时间序列却又是十分相似的。例如正弦函数  和余弦函数 ，只是平移了  个长度而已。本文将会介绍一些基于形状的时间序列的距离算法，并且介绍如何在给定时间序列的情况下，在时间序列数据库中寻找相似的时间序列。

基于动态规划的相似度计算方法

DTW 的经典算法

基于动态规划的相似度计算的典型代表就是 DTW（Dynamical Time Warping ）和 Frechet 距离。在这里将会主要介绍 DTW 算法。详细来说，DTW 算法是为了计算语音信号处理中，由于两个人说话的时长不一样，但是却是类似的一段话，欧几里德算法不完全能够解决这类问题。在这种情况下，DTW 算法就被发展出来。DTW 算法是为了计算两条时间序列的最佳匹配点，假设我们有两条时间序列 Q 和 C，长度都是 n，并且  和 。首先我们可以建立一个  的矩阵， 位置的元素是 ，这里的 dist 可以使用  范数。其次，我们想找到一条路径，使得这个矩阵的累积距离最小，而这条路则是两条时间序列之间的最佳匹配。在这里，我们可以假设这条路径是 ，其中  的每个元素表示时间序列 Q 中的第 i 个元素和时间序列 C 中的第 j 个元素之间的距离. i.e. 。

现在我们需要找到一条路径使得

.

这条路径就是动态规划的解，它满足一个动态规划方程：对于 ，有

其初始状态是   最终的取值  就是我们需要的解，也就是两条时间序列的 DTW 距离。按照上面的算法，DTW 算法的时间复杂度是 。特别地，

1. 如果  时，则  表示最后的距离；
2. 如果  时，则  表示最后的距离。
3. 如果  时，则  表示最后的距离。

Remark. 

DTW 算法不满足三角不等式。例如：，则

DTW 的加速算法

某些时候，我们可以添加一个窗口长度的限制，换言之，如果要比较  与  的话，i 与 j 需要满足 ，这里的 w 表示窗口长度。因此算法的描述如下：

初始条件和之前一样。

这里的  取值范围是：对每一个 ，需要  满足。

相似时间序列的搜索

相似的时间序列的搜索问题一般是这样的：

Question 1. 给定一条时间序列  和一个时间序列的数据库 。通过某种相似度或者距离计算方法，计算出给定的时间序列  与时间序列数据库中  中最相似的时间序列。

Question 2. 给定一条时间序列  和一个时间序列的数据库 ，以及正整数 。从数据库  中寻找与给定的时间序列  最相似的  条时间序列。

DTW 算法的下界 LB_Kim

对于两条时间序列 Q 和 C 而言，可以分别提取它们的四个特征，那就是最大值，最小值，第一个元素的值，最后一个元素的值。这样可以计算出 LB_Kim

可以证明 .

DTW 算法的下界 LB_Yi

有学者在 LB_Kim 的基础上提出了一种下界的计算方法，那就是 LB_Yi，有兴趣的读者可以参见：efficient retrieval of similar time sequences under time warping， 1998.

DTW 算法的下界 LB_Keogh

但是，如果是基于 DTW 的距离计算方法，每次都要计算两条时间序列的 DTW 距离，时间复杂度是 。于是就有学者研究是否存在 DTW 距离的下界表示，也就是说找到一个合适的下界，Lower Bound of DTW。每次判断 Lower Bound of DTW 是否小于当前的最小距离，如果下界高于最小距离，就不需要进行 DTW 的计算；否则开始计算 DTW 的值。如果下界的计算速度足够快，并且下界足够精准的话，可以大量的压缩搜索的时间。于是，Keogh 提出了下界的计算方法。

（1）首先定义时间序列 Q 的上下界。i.e. ，给定一个窗口的取值 r，得到 ，。

（2）定义公式：

其中，LBKeogh 满足性质：

对于任意两条长度为 n 的时间序列 Q 和 C，对于任意的窗口 ，有不等式  成立。

所以，可以把搜索的算法进行加速：

Lower Bounding Sequential Scan(Q)
    best_so_far = infinity 
    for all sequences in database
        LB_dist = lower_bound_distance(C_{i},Q)
        if LB_dist < best_so_far
            true_dist = DTW(C_{i},Q)
            if true_dist < best_so_far
                best_so_far = true_dist
                index_of_best_match = i
            endif
        endif
    endfor

在论文 Exact Indexing of Dynamic Time Warping 里面，作者还尝试将 Piecewise Constant Approximation 与 LB_Keogh 结合起来，提出了 LB_PAA 的下界。它满足条件 .

总结

本文初步介绍了 DTW 算法以及它的下界算法，包括 LB_Kim, LB_Keogh 等，以及时间序列数据库的搜索算法。

Riemann 积分和 Lebesgue 积分是数学中两个非常重要的概念。本文将会从 Riemann 积分和 Lebesgue 积分的定义出发，介绍它们各自的性质和联系。

积分

Riemann 积分

Riemann 积分虽然被称为 Riemann 积分，但是在 Riemann 之前就有学者对这类积分进行了详细的研究。早在阿基米德时代，阿基米德为了计算曲线  在 [0,1] 区间上与 X 坐标轴所夹的图形面积，就使用了 Riemann 积分的思想。 他把 [0,1] 区间等长地切割成 n 段，每一段使用一个长方形去逼近  这条曲线的分段面积，再把 n 取得很大，所以得到当 n 趋近于无穷的时候，就知道该面积其实是 1/3。

下面来看一下 Riemann 积分的详细定义。

考虑定义在闭区间 [a,b] 上的函数 ，

取一个有限的点列 ， 表示这些区间长度的最大值，在这里 。在每一个子区间上 上取出一个点 。而函数  关于以上取样分割的 Riemann 和就是以下公式：

当我们说该函数  在闭区间  上的取值是  的意思是：

对于任意 ，存在  使得对于任意取样分割，当  时，就有

通常来说，用符号来表示就是：.

用几幅图来描述 Riemann 积分的思想就是：

Lebesgue 积分

Riemann 积分是为了计算曲线与 X 轴所围成的面积，而 Lebesgue 积分也是做同样的事情，但是计算面积的方法略有不同。要想直观的解释两种积分的原理，可以参见下图：

Riemann 积分是把一条曲线的底部分成等长的区间，测量每一个区间上的曲线高度，所以总面积就是这些区间与高度所围成的面积和。

Lebesgue 积分是把曲线化成等高线图，每两根相邻等高线的差值是一样的。每根等高线之内含有它所圈着的长度，因此总面积就是这些等高线内的面积之和。

用再形象一点的语言来描述就是：吃一块汉堡有多种方式

1. Riemann 积分：从一个角落开始一口一口吃，每口都包含所有的配料；
2. Lebesgue 积分：从最上层开始吃，按照“面包-配菜-肉-蛋-面包”的节奏，一层一层来吃。

再看一幅图的表示就是：Riemann 积分是按照蓝色的数字顺序相加的，Lebesgue 积分是按照红色的数字来顺序相加的。

基于这些基本的思想，就可以给出 Lebesgue 积分的定义：

简单函数指的是对指示函数的有限线性组合，i.e. ，这里的  是系数， 是可测集合， 表示指示函数。当  非负时，令

如果  是一个非负可测函数时，可以定义函数  在可测集合  上的 Lebesgue 积分是：

,

这里的  指的是非负简单函数， 表示零函数，这里的大小关系表示对定义域内的每个点都要成立。

而对于可测函数时，可以把可测函数  转换成 ，而这里的  和  都是非负可测函数。所以可以定义任意可测函数的 Lebesgue 积分如下：

.

Riemann 积分与Lebesgue 积分的关系

定义了两种积分之后，也许有人会问它们之间是否存在矛盾？其实，它们之间是不矛盾的，因为有学者证明了这样的定理：

如果有界函数  在闭区间  是 Riemann 可积的，则它也是 Lebesgue 可积的，并且它们的积分值相等：

.

左侧是表示 Riemann 积分，右侧表示 Lebesgue 积分。

用形象化一点的语言描述就是：无论从角落一口一口地吃汉堡，还是从顶至下一层一层吃，所吃的汉堡都是同一个。

但是 Lebesgue 积分比 Riemann 积分有着更大的优势，例如 Dirichlet 函数，

1. 当  是有理数时，；
2. 当  是无理数时，.

Dirichlet 函数是定义在实数轴的函数，并且值域是 ，无法画出函数图像，它不是 Riemann 可积的，但是它 Lebesgue 可积。

时间序列

提到时间序列，也就是把以上所讨论的连续函数换成离散函数而已，把定义域从一个闭区间  换成  这样的定义域而已。所以，之前所讨论的很多连续函数的想法都可以应用在时间序列上。

时间序列的表示 — 基于 Riemann 积分

现在我们可以按照 Riemann 积分的计算方法来表示一个时间序列的特征，于是就有学者把时间序列按照横轴切分成很多段，每一段使用某个简单函数（线性函数等）来表示，于是就有了以下的方法：

1. 分段线性逼近（Piecewise Linear Approximation）
2. 分段聚合逼近（Piecewise Aggregate Approximation）
3. 分段常数逼近（Piecewise Constant Approximation）

说到这几种算法，其实最本质的思想就是进行数据降维的工作，用少数的数据来进行原始时间序列的表示（Representation）。用数学化的语言来描述时间序列的数据降维（Data Reduction）就是：把原始的时间序列  用  来表示，其中 。那么后者就是原始序列的一种表示（representation）。

分段聚合逼近（Piecewise Aggregate Approximation）— 类似 Riemann 积分

在这种算法中，分段聚合逼近（Piecewise Aggregate Approximation）是一种非常经典的算法。假设原始的时间序列是 ，定义 PAA 的序列是：，

其中

.

在这里 。用图像来表示那就是：

至于分段线性逼近（Piecewise Linear Approximation）和分段常数逼近（Piecewise Constant Approximation），只需要在  的定义上稍作修改即可。

符号特征（Symbolic Approximation）— 类似用简单函数来计算 Lebesgue 积分

在推荐系统的特征工程里面，特征通常来说可以做归一化，二值化，离散化等操作。例如，用户的年龄特征，一般不会直接使用具体的年月日，而是划分为某个区间段，例如 0~6（婴幼儿时期），7~12（小学），13~17（中学），18~22（大学）等阶段。

其实在得到分段特征之后，分段特征在某种程度上来说依旧是某些连续值，能否把连续值划分为一些离散的值呢？于是就有学者使用一些符号来表示时间序列的关键特征，也就是所谓的符号表示法（Symbolic Representation）。下面来介绍经典的 SAX Representation。

如果我们希望使用  个符号来表示时间序列，那么我们其实可以考虑正态分布 ，用 来表示 Gauss 曲线下方的一些点，而这些点把 Gauss 曲线下方的面积等分成了  段。用  表示  个字母。

SAX 方法的流程如下：

1. 正规化（normalization）：把原始的时间序列映射到一个新的时间序列，新的时间序列满足均值为零，方差为一的条件。
2. 分段表示（PAA）：。
3. 符号表示（SAX）：如果 ，那么 ；如果 ，那么 ，在这里 ；如果 ，那么 。

于是，我们就可以用  这  个字母来表示原始的时间序列了。

时间序列的表示 — 基于 Lebesgue 积分

要想考虑一个时间序列的值分布情况，其实就类似于 Lebesgue 积分的计算方法，考虑它们的分布情况，然后使用某些函数去逼近时间序列。要考虑时间序列的值分布情况，可以考虑熵的概念。

熵（Entropy）

通常来说，要想描述一种确定性与不确定性，熵（entropy）是一种不错的指标。对于离散空间而言，一个系统的熵（entropy）可以这样来表示：

.

如果一个系统的熵（entropy）越大，说明这个系统就越混乱；如果一个系统的熵越小，那么说明这个系统就更加确定。

提到时间序列的熵特征，一般来说有几个经典的熵指标，其中有一个就是 binned entropy。

分桶熵（Binned Entropy）

从熵的定义出发，可以考虑把时间序列的值进行分桶的操作，例如，可以把 [min, max] 这个区间等分为十个小区间，那么时间序列的取值就会分散在这十个桶中。根据这个等距分桶的情况，就可以计算出这个概率分布的熵（entropy）。i.e. Binned Entropy 就可以定义为：

其中  表示时间序列 X 的取值落在第 k 个桶的比例（概率），maxbin 表示桶的个数，len(X) 表示时间序列 X 的长度。

如果一个时间序列的 Binned Entropy 较大，说明这一段时间序列的取值是较为均匀的分布在 [min, max] 之间的；如果一个时间序列的 Binned Entropy 较小，说明这一段时间序列的取值是集中在某一段上的。

总结

在本篇文章中，笔者从 Riemann 积分和 Lebesgue 积分出发，介绍了它们的基本概念，性质和联系。然后从两种积分出发，探讨了时间序列的分段特征，时间序列的熵特征。在未来的 Blog 中，笔者将会介绍时间序列的更多相关内容。

在文本挖掘中，可以通过 Word2Vec 生成的向量以及向量的内积，或者根据语义和词性来判断两个词语是否是近义词。在时间序列的挖掘中，同样可以找到一些方法来描述两条时间序列是否相似。

在介绍时间序列的距离之前，笔者感觉需要回顾一下数学中度量空间和内积空间的定义。

度量空间

在数学里面，集合  上的距离函数定义为 ，其中  表示实数集合，并且函数  满足以下几个条件：

1. ，并且  当且仅当 ;
2. ，也就是满足对称性；
3. ，也就是三角不等式。

满足这三个条件的距离函数称为度量，具有某种度量的集合则叫做度量空间。

Remark.

在本文下面的时间序列距离的各种定义中，这些距离不一定满足三角不等式。例如动态时间算法（DTW）就不满足三角不等式。

内积空间

一个内积空间是域 （其中  或者 ）上的向量空间 与一个内积（映射）所构成， 上的一个内积定义为正定，非退化的共轭双线性形式，记成 ，它满足以下设定：

1. 对于任意的 ，有 

2. 共轭双线性形式指的是：

3. 非负性：

4. 非退化：从  到对偶空间  的映射： 是同构映射。

基于欧几里德距离的相似度计算

假设两条时间序列曲线为  和 ，于是可以使用欧几里德空间里面的  范数来表示两个时间序列之间的距离。用公式来描述就是：

基于相关性的相似度计算方法

Pearson 系数（Pearson Coefficient）

其中，

, 。

Pearson 系数的性质如下：

1. 如果两条时间序列 ，则  表是它们是完全一致的，如果两条时间序列 ，则  表示它们之间是负相关的。
2. .

可以基于 Pearson 系数来制定两条时间序列之间的距离：

其中.

The First Order Temporal Correlation Coefficient

这个相关性系数与 Pearson 系数类似，但是略有不同，其定义为：

 的性质：

1. 
2.  表示两条时间序列持有类似的趋势， 它们会同时上涨或者下跌，并且涨幅或者跌幅也是类似的。
3.  表示两条时间序列的上涨和下跌趋势恰好相反。
4.  表示两条时间序列在单调性方面没有相关性。

基于 CORT，同样可以定义时间序列的距离，用公式描述如下：

其中， 可以用  来计算，而

是一个递减函数。

基于自相关系数的距离（Autocorrelation-based distance)

假设时间序列是 ，对于任意的 ，可以定义自相关系数为：

,

其中  分别表示该时间序列的均值和方差。该公式相当于是比较整个时间序列  的两个子序列的相似度（Pearson 系数），这两个子序列分别是  和 。

于是，通过给定一个正整数 ，可以对每一个时间序列得到一组自相关系数的向量，用公式描述如下：

对于  的情况，可以假定  和 。于是，可以定义时间序列之间的距离如下：

.

其中的  表示一个  的矩阵。它有着很多种选择，例如：

（1） 表示单位矩阵。用公式表示就是

.

（2） 表示一个 的对角矩阵，其中 。此时相当于一个带权重的求和公式。

.

除了自相关系数（Autocorrelation Coefficients）之外，也可以考虑偏自相关系数（Partial Autocorrelation Coefficients），使用 PACFs 来取代 ACFs。这样，使用同样的定义方式就可以得到  和  两个距离公式。

基于周期性的相似度计算方法

这里会介绍基于周期图表（Periodogram-based）的距离计算方法。其大体思想就是通过 Fourier 变换得到一组参数，然后通过这组参数来反映原始的两个时间序列时间的距离。用数学公式来描述就是：

,

.

其中 ，，。这里的  表示 Gauss 取整函数。

（1）用原始的特征来表示距离：

.

（2）用正则化之后的特征来描述就是：

,

其中 ，，和  表示  的标准差（sample variance）。

（3）用取对数之后的特征表示：

.

基于模型的相似度计算

Piccolo 距离

基于模型的相似度判断本质上是用一个模型和相应的一组参数去拟合某条时间序列，然后得到最优的一组参数，计算两个时间序列所得到的最优参数的欧几里德距离即可。

 模型有自己的 AR 表示，因此可以得到相应的一组参数 ，所以，对于每一条时间序列，都可以用一组最优的参数去逼近。如果

分别表示  和  对于时间序列  和  的参数估计，则 Piccolo 距离如下：

,

其中 ， 当 ，并且  当 。 当 ，并且  当 。

Maharaj 距离

按照之前的描述，可以增加一个矩阵来修改 Piccolo 距离：

其中  和  表示  模型对于  和  的参数估计，和 Piccolo 距离一样。， 和  表示时间序列的方差， 和  表示时间序列的 sample covariance 矩阵。

基于 Cepstral 的距离

考虑时间序列  满足  的结构，i.e. ，这里的  表示 AR 模型的参数， 表示白噪声（均值为 0，方差为 1 的 Gauss 正态分布）。于是可以从这些参数定义 LPC 系数如下：

，

 当 ，

 当 。

所以，LPC 的距离定义是：

.

总结

在本文中，介绍了时间序列之间距离的计算方法，包括基于  范数的距离，基于相关性的距离，基于周期图表的计算方法，基于模型的计算方法。

提到时间序列，大家能够想到的就是一串按时间排序的数据，但是在这串数字背后有着它特殊的含义，那么如何进行时间序列的表示（Representation），如何进行时间序列的信息提取（Information Extraction）就成为了时间序列研究的关键问题。

就笔者的个人经验而言，其实时间序列的一些想法和文本挖掘是非常类似的。通常来说句子都是由各种各样的词语组成的，并且一般情况下都是“主谓宾”的句子结构。于是就有人希望把词语用一个数学上的向量描述出来，那么最经典的做法就是使用 one – hot 的编码格式。i.e. 也就是对字典里面的每一个词语进行编码，一个词语对应着一个唯一的数字，例如 0，1，2 这种形式。one hot 的编码格式是这行向量的长度是词典中词语的个数，只有一个值是1，其余的取值是0，也就是 (0,…,0,1,0,…,0) 这种样子。但是在一般情况下，词语的个数都是非常多的，如何使用一个维度较小的向量来表示一个词语就成为了一个关键的问题。几年前，GOOGLE 公司开源了 Word2vec 开源框架，把每一个词语用一串向量来进行描述，向量的长度可以自行调整，大约是100~1000 不等，就把原始的 one-hot 编码转换为了一个低维空间的向量。在这种情况下，机器学习的很多经典算法，包括分类，回归，聚类等都可以在文本上得到巨大的使用。Word2vec 是采用神经网络的思想来提取每个词语与周边词语的关系，从而把每个词语用一个低维向量来表示。在这里，时间序列的特征提取方法与 word2vec 略有不同，后面会一一展示这些技巧。

时间序列的统计特征

提到时间序列的统计特征，一般都能够想到最大值（max），最小值（min），均值（mean），中位数（median），方差（variance），标准差（standard variance）等指标，不过一般的统计书上还会介绍两个指标，那就是偏度（skewness）和峰度（kuriosis）。如果使用时间序列  来表示长度为  的时间序列，那么这些统计特征用数学公式来表示就是：

其中  和  分别表示时间序列  的均值和方差。

时间序列的熵特征

为什么要研究时间序列的熵呢？请看下面两个时间序列：

时间序列（1）：(1,2,1,2,1,2,1,2,1,2,…)

时间序列（2）：(1,1,2,1,2,2,2,2,1,1,…)

在时间序列（1）中，1 和 2 是交替出现的，而在时间序列（2）中，1 和 2 是随机出现的。在这种情况下，时间序列（1）则更加确定，时间序列（2）则更加随机。并且在这种情况下，两个时间序列的统计特征，例如均值，方差，中位数等等则是几乎一致的，说明用之前的统计特征并不足以精准的区分这两种时间序列。

通常来说，要想描述一种确定性与不确定性，熵（entropy）是一种不错的指标。对于离散空间而言，一个系统的熵（entropy）可以这样来表示：

如果一个系统的熵（entropy）越大，说明这个系统就越混乱；如果一个系统的熵越小，那么说明这个系统就更加确定。

提到时间序列的熵特征，一般来说有几个经典的例子，那就是 binned entropy，approximate entropy，sample entropy。下面来一一介绍时间序列中这几个经典的熵。

Binned Entropy

从熵的定义出发，可以考虑把时间序列  的取值进行分桶的操作。例如，可以把  这个区间等分为十个小区间，那么时间序列的取值就会分散在这十个桶中。根据这个等距分桶的情况，就可以计算出这个概率分布的熵（entropy）。i.e. Binned Entropy 就可以定义为：

其中  表示时间序列  的取值落在第  个桶的比例（概率）， 表示桶的个数， 表示时间序列  的长度。

如果一个时间序列的 Binned Entropy 较大，说明这一段时间序列的取值是较为均匀的分布在  之间的；如果一个时间序列的 Binned Entropy 较小，说明这一段时间序列的取值是集中在某一段上的。

Approximate Entropy

回到本节的问题，如何判断一个时间序列是否具备某种趋势还是随机出现呢？这就需要介绍 Approximate Entropy 的概念了，Approximate Entropy 的思想就是把一维空间的时间序列提升到高维空间中，通过高维空间的向量之间的距离或者相似度的判断，来推导出一维空间的时间序列是否存在某种趋势或者确定性。那么，我们现在可以假设时间序列  的长度是 ，同时 Approximate Entropy 函数拥有两个参数， 与 ，下面来详细介绍 Approximate Entropy 的算法细节。

Step 1. 固定两个参数，正整数  和正数 ，正整数  是为了把时间序列进行一个片段的提取，正数  是表示时间序列距离的某个参数。i.e. 需要构造新的  维向量如下：

Step 2. 通过新的向量 ，可以计算出哪些向量与  较为相似。i.e.

在这里，距离  可以选择  范数。在这个场景下，距离  通常选择为  范数。

Step 3. 考虑函数

Step 4. Approximate Entropy 可以定义为：

Remark.

1. 正整数  一般可以取值为 2 或者 3， 会基于具体的时间序列具体调整；
2. 如果某条时间序列具有很多重复的片段（repetitive pattern）或者自相似性（self-similarity pattern），那么它的 Approximate Entropy 就会相对小；反之，如果某条时间序列几乎是随机出现的，那么它的 Approximate Entropy 就会相对较大。

Sample Entropy

除了 Approximate Entropy，还有另外一个熵的指标可以衡量时间序列，那就是 Sample Entropy，通过自然对数的计算来表示时间序列是否具备某种自相似性。

按照以上 Approximate Entropy 的定义，可以基于  与  定义两个指标  和 ，分别是

其中， 表示集合的元素个数。根据度量 （无论是 ）的定义可以知道，因此 Sample Entropy 总是非负数，i.e.

Remark.

1. Sample Entropy 总是非负数；
2. Sample Entropy 越小表示该时间序列具有越强的自相似性（self similarity）。
3. 通常来说，在 Sample Entropy 的参数选择中，可以选择 .

时间序列的分段特征

即使时间序列有一定的自相似性（self-similarity），能否说明这两条时间序列就完全相似呢？其实答案是否定的，例如：两个长度都是 1000 的时间序列，

时间序列（1）： [1,2] * 500

时间序列（2）： [1,2,3,4,5,6,7,8,9,10] * 100

其中，时间序列（1）是 1 和 2 循环的，时间序列（2）是 1~10 这样循环的，它们从图像上看完全是不一样的曲线，并且它们的 Approximate Entropy  和 Sample Entropy 都是非常小的。那么问题来了，有没有办法提炼出信息，从而表示它们的不同点呢？答案是肯定的。

首先，我们可以回顾一下 Riemann 积分和 Lebesgue 积分的定义和不同之处。按照下面两幅图所示，Riemann 积分是为了算曲线下面所围成的面积，因此把横轴划分成一个又一个的小区间，按照长方形累加的算法来计算面积。而 Lebesgue 积分的算法恰好相反，它是把纵轴切分成一个又一个的小区间，然后也是按照长方形累加的算法来计算面积。

之前的 Binned Entropy 方案是根据值域来进行切分的，好比 Lebesgue 积分的计算方法。现在我们可以按照 Riemann 积分的计算方法来表示一个时间序列的特征，于是就有学者把时间序列按照横轴切分成很多段，每一段使用某个简单函数（线性函数等）来表示，于是就有了以下的方法：

1. 分段线性逼近（Piecewise Linear Approximation）
2. 分段聚合逼近（Piecewise Aggregate Approximation）
3. 分段常数逼近（Piecewise Constant Approximation）

说到这几种算法，其实最本质的思想就是进行数据降维的工作，用少数的数据来进行原始时间序列的表示（Representation）。用数学化的语言来描述时间序列的数据降维（Data Reduction）就是：把原始的时间序列  用  来表示，其中 。那么后者就是原始序列的一种表示（representation）。

分段聚合逼近（Piecewise Aggregate Approximation）— 类似 Riemann 积分

在这种算法中，分段聚合逼近（Piecewise Aggregate Approximation）是一种非常经典的算法。假设原始的时间序列是 ，定义 PAA 的序列是：，

其中

.

在这里 。用图像来表示那就是：

至于分段线性逼近（Piecewise Linear Approximation）和分段常数逼近（Piecewise Constant Approximation），只需要在  的定义上稍作修改即可。

符号逼近（Symbolic Approximation）— 类似 Riemann 积分

在推荐系统的特征工程里面，特征通常来说可以做归一化，二值化，离散化等操作。例如，用户的年龄特征，一般不会直接使用具体的年月日，而是划分为某个区间段，例如 0~6（婴幼儿时期），7~12（小学），13~17（中学），18~22（大学）等阶段。

其实在得到分段特征之后，分段特征在某种程度上来说依旧是某些连续值，能否把连续值划分为一些离散的值呢？于是就有学者使用一些符号来表示时间序列的关键特征，也就是所谓的符号表示法（Symbolic Representation）。下面来介绍经典的 SAX Representation。

如果我们希望使用  个符号，例如  来表示时间序列。同时考虑正态分布 ，用 来表示 Gauss 曲线下方的一些点，而这些点把 Gauss 曲线下方的面积等分成了  段。

SAX 方法的流程如下：

Step 1. 正规化（normalization）：也就是该时间序列被映射到均值为零，方差为一的区间内。

Step 2. 分段表示（PAA）：。

Step 3. 符号表示（SAX）：如果 ，那么 ；如果 ，那么 ；如果 ，那么 。

于是，我们就可以用  这  个字母来表示原始的时间序列了。

、

总结

在本篇文章中，我们介绍了时间序列的一些表示方法（Representation），其中包括时间序列统计特征，时间序列的熵特征，时间序列的分段特征。在下一篇文章中，我们将会介绍时间序列的相似度计算方法。

https://machinelearningmastery.com/convert-time-series-supervised-learning-problem-python/

Machine learning methods like deep learning can be used for time series forecasting.

Before machine learning can be used, time series forecasting problems must be re-framed as supervised learning problems. From a sequence to pairs of input and output sequences.

In this tutorial, you will discover how to transform univariate and multivariate time series forecasting problems into supervised learning problems for use with machine learning algorithms.

After completing this tutorial, you will know:

• How to develop a function to transform a time series dataset into a supervised learning dataset.
• How to transform univariate time series data for machine learning.
• How to transform multivariate time series data for machine learning.

Let’s get started.

Time Series vs Supervised Learning

Before we get started, let’s take a moment to better understand the form of time series and supervised learning data.

A time series is a sequence of numbers that are ordered by a time index. This can be thought of as a list or column of ordered values.

For example:

A supervised learning problem is comprised of input patterns (X) and output patterns (y), such that an algorithm can learn how to predict the output patterns from the input patterns.

For example:

For more on this topic, see the post:

• Time Series Forecasting as Supervised Learning

Pandas shift() Function

A key function to help transform time series data into a supervised learning problem is the Pandas shift() function.

Given a DataFrame, the shift() function can be used to create copies of columns that are pushed forward (rows of NaN values added to the front) or pulled back (rows of NaN values added to the end).

This is the behavior required to create columns of lag observations as well as columns of forecast observations for a time series dataset in a supervised learning format.

Let’s look at some examples of the shift function in action.

We can define a mock time series dataset as a sequence of 10 numbers, in this case a single column in a DataFrame as follows:

Running the example prints the time series data with the row indices for each observation.

We can shift all the observations down by one time step by inserting one new row at the top. Because the new row has no data, we can use NaN to represent “no data”.

The shift function can do this for us and we can insert this shifted column next to our original series.

Running the example gives us two columns in the dataset. The first with the original observations and a new shifted column.

We can see that shifting the series forward one time step gives us a primitive supervised learning problem, although with X and y in the wrong order. Ignore the column of row labels. The first row would have to be discarded because of the NaN value. The second row shows the input value of 0.0 in the second column (input or X) and the value of 1 in the first column (output or y).

We can see that if we can repeat this process with shifts of 2, 3, and more, how we could create long input sequences (X) that can be used to forecast an output value (y).

The shift operator can also accept a negative integer value. This has the effect of pulling the observations up by inserting new rows at the end. Below is an example:

Running the example shows a new column with a NaN value as the last value.

We can see that the forecast column can be taken as an input (X) and the second as an output value (y). That is the input value of 0 can be used to forecast the output value of 1.

Technically, in time series forecasting terminology the current time (t) and future times (t+1, t+n) are forecast times and past observations (t-1, t-n) are used to make forecasts.

We can see how positive and negative shifts can be used to create a new DataFrame from a time series with sequences of input and output patterns for a supervised learning problem.

This permits not only classical X -> y prediction, but also X -> Y where both input and output can be sequences.

Further, the shift function also works on so-called multivariate time series problems. That is where instead of having one set of observations for a time series, we have multiple (e.g. temperature and pressure). All variates in the time series can be shifted forward or backward to create multivariate input and output sequences. We will explore this more later in the tutorial.

The series_to_supervised() Function

We can use the shift() function in Pandas to automatically create new framings of time series problems given the desired length of input and output sequences.

This would be a useful tool as it would allow us to explore different framings of a time series problem with machine learning algorithms to see which might result in better performing models.

In this section, we will define a new Python function named series_to_supervised() that takes a univariate or multivariate time series and frames it as a supervised learning dataset.

The function takes four arguments:

• data: Sequence of observations as a list or 2D NumPy array. Required.
• n_in: Number of lag observations as input (X). Values may be between [1..len(data)] Optional. Defaults to 1.
• n_out: Number of observations as output (y). Values may be between [0..len(data)-1]. Optional. Defaults to 1.
• dropnan: Boolean whether or not to drop rows with NaN values. Optional. Defaults to True.

The function returns a single value:

• return: Pandas DataFrame of series framed for supervised learning.

The new dataset is constructed as a DataFrame, with each column suitably named both by variable number and time step. This allows you to design a variety of different time step sequence type forecasting problems from a given univariate or multivariate time series.

Once the DataFrame is returned, you can decide how to split the rows of the returned DataFrame into X and y components for supervised learning any way you wish.

The function is defined with default parameters so that if you call it with just your data, it will construct a DataFrame with t-1 as X and t as y.

The function is confirmed to be compatible with Python 2 and Python 3.

The complete function is listed below, including function comments.

Can you see obvious ways to make the function more robust or more readable?
Please let me know in the comments below.

Now that we have the whole function, we can explore how it may be used.

One-Step Univariate Forecasting

It is standard practice in time series forecasting to use lagged observations (e.g. t-1) as input variables to forecast the current time step (t).

This is called one-step forecasting.

The example below demonstrates a one lag time step (t-1) to predict the current time step (t).

Running the example prints the output of the reframed time series.

We can see that the observations are named “var1” and that the input observation is suitably named (t-1) and the output time step is named (t).

We can also see that rows with NaN values have been automatically removed from the DataFrame.

We can repeat this example with an arbitrary number length input sequence, such as 3. This can be done by specifying the length of the input sequence as an argument; for example:

The complete example is listed below.

Again, running the example prints the reframed series. We can see that the input sequence is in the correct left-to-right order with the output variable to be predicted on the far right.

Multi-Step or Sequence Forecasting

A different type of forecasting problem is using past observations to forecast a sequence of future observations.

This may be called sequence forecasting or multi-step forecasting.

We can frame a time series for sequence forecasting by specifying another argument. For example, we could frame a forecast problem with an input sequence of 2 past observations to forecast 2 future observations as follows:

The complete example is listed below:

Running the example shows the differentiation of input (t-n) and output (t+n) variables with the current observation (t) considered an output.

Multivariate Forecasting

Another important type of time series is called multivariate time series.

This is where we may have observations of multiple different measures and an interest in forecasting one or more of them.

For example, we may have two sets of time series observations obs1 and obs2 and we wish to forecast one or both of these.

We can call series_to_supervised() in exactly the same way.

For example:

Running the example prints the new framing of the data, showing an input pattern with one time step for both variables and an output pattern of one time step for both variables.

Again, depending on the specifics of the problem, the division of columns into X and Y components can be chosen arbitrarily, such as if the current observation of var1 was also provided as input and only var2was to be predicted.

You can see how this may be easily used for sequence forecasting with multivariate time series by specifying the length of the input and output sequences as above.

For example, below is an example of a reframing with 1 time step as input and 2 time steps as forecast sequence.

Running the example shows the large reframed DataFrame.

Experiment with your own dataset and try multiple different framings to see what works best.

Summary

In this tutorial, you discovered how to reframe time series datasets as supervised learning problems with Python.

Specifically, you learned:

• About the Pandas shift() function and how it can be used to automatically define supervised learning datasets from time series data.
• How to reframe a univariate time series into one-step and multi-step supervised learning problems.
• How to reframe multivariate time series into one-step and multi-step supervised learning problems.

Do you have any questions?
Ask your questions in the comments below and I will do my best to answer.

论文：Exact Discovery of Time Series Motifs

Speeded up Brute Force Motif Discovery:

Github：https://github.com/saifuddin778/mkalgo

但是感觉有一行比较奇怪，应该是 ，而不是，因为  是递增排列的，并且 best-so-far > 0.

Generalization to multiple reference points:

https://github.com/nicholasg3/motif-mining/tree/95bbb05ac5d0f9e90134a67a789ea7e607f22cea

注意：

for j = 1 to m-offset 而不是 for j = 1 to R

Time Series Clustering with Dynamic Time Warping (DTW)

https://github.com/goodmattg/wikipedia_kaggle