时间序列相似性属于曲线相似性/曲线匹配(curve matching)领域的内容，在这一领域，有许多有用的方法，但是国内的博客上鲜有这方面的内容，因此我选取了几种常用的方法进行一下综述性的阐述。

衡量相似性之前，我们首先定义“相似”。

正常情况下，我们认为x,y,z是形状相似的，在这三条曲线中，我们认为y，z是最相似的两条曲线(因为y,z的距离最近)。

ok，那我们先来看看寻常意义上的相似：距离最近且形状相似。本文主要详细介绍时间序列相似度计算的DTW算法和PLR算法。

1. 欧式距离

要衡量距离与形状，显然欧式距离是一个天然完美的指标，上图中我们正是基于欧式距离认为y与z是最相似的，欧式距离在诸多算法都有广泛的应用。对于长度相同的序列，计算每两点之间的距离然后求和，距离越小相似度越高(whole matching)。对于不同长度的序列，一般有两种方法处理：

1)子序列匹配(subsequence matching): 找出长序列中与短序列最相似的部分。举个栗子，设序列

,序列

，其中

。滚动地计算A与B的距离：

，然后找出所有d中的最小值，该距离所对应的A序列的索引即为A中与B最相似的部分。

2)滑动窗口：微软在2001年在Dimensionality Reduction for Fast Similarity Search文中提出为了减少算法复杂度，可以复制B序直到与A序列等长。

由于微软之后使用了独特的降维方法，且计算复杂度不是本文考虑的主要内容，因此，在涉及长短序列相似度计算的时候，本文均使用第一种方法。

似乎时间序列的相似性度量的计算可以就此为止了，然而远非如此。

天津大学的XIAO-LI DONG, CHENG-KUI GU, ZHENG-OU WANG在2006年Research on shape-based time series similarity measure[C]//2006 International Conference on Machine Learning and Cybernetics. IEEE, 2006: 1253-1258一文中指出了欧式距离用于衡量时间序列相似性的三个缺陷：不能辨别形状相似性

不能反映趋势动态变化幅度的相似性

基于点距离的计算不能反映不同分析频率的不同

举个栗子：

A与B的变化趋势几乎完全相反，A与C的变化趋势几乎完全相同。如果使用欧式距离去度量，那么结论就是A与B是最相似的。而实际上，在变化是A与C是相似的。

为了进一步加强对欧式距离的理解，我们不妨再举一个简单的例子：

正常来说，我们认为与y1最相似的是y3，实际上，y3就是y1向下平移得到的。然而欧式距离告诉我们，距离y1最近的是y2。

下面是使用Python进行模拟的源代码：

import numpy as np

import matplotlib.pyplot as plt

x=np.arange(0,np.pi*2,0.1)

y1=np.sin(x)

y2=np.cos(x)-2

y3=y1-2

plt.plot(y1)

plt.plot(y2)

plt.plot(y3)

plt.legend(['y1','y2','y3'])

def dis(x,y):

return(sum(x-y)**2)

dis(y1,y2)

Out[15]: 15831.914509398732

dis(y1,y3)

Out[16]: 15876.0

欧式距离对形状的度量如此糟糕，有没有更好的模型能度量形状呢？

2. 模式距离Pattern distance

首先引入一个算法，PLR(piecewise linear representation)分段线性表示。

一个时间序列，无非有三种状态：上升、下降和不变，我们将这种状态对应表示为

。假设有某个长度为S的序列，我们将其等分为K段。对于每一段计算一个斜率，斜率为正表示上升，为负表示下降，为0表示不变。

我画了一个草图，可能emm……不是很美观。

那么我们就可以将这个序列表示为[1,1,0,-1...]这样的序列，将相邻的相同模式进行合并，我们得到[1,0,-1...]的序列。这就是PLR算法。

关于PLR算法的点的分割，为了便于说理，我这里直接用的是等分的方式，然而观察上图可知，第三个模式表示为了0，实际上第三个模式是一个尖峰，是先上升后下降的，所以等分切割的方法并不科学。

KEOGH E. 在Fast similarity search in the presence of longitudinal scaling in time series data bases中提出的自底向上的搜索方法很好的解决了这个问题，感兴趣的读者可以自行了解。

因为我们将相邻的相同模式进行了合并，所以我们得到的模式序列一定是1,-1,0间隔排列的，每一个模式可能跨越了不同的时间长度，在合并模式后，序列S1可能有N个模式，S2可能有M个模式。现在我们要将他们等模式数化。

如上图所示，经过PLR后，我们将S1，S2表示为：

设

分表表示S1,S2的第i,j个模式，即

，t表示时间，且无论怎么切割，最后的终点都是相等的，即

。结合上图，具体一点就是：

等模式数化后将他们变形为：

说白了就让他们使用共同的分割点，以获得最后长度相等的模式序列。

OK现在我们已经完全定义好了何为模式，接下来我们计算距离。这里我们使用绝对值距离，当然使用欧式距离也是可以的。模式距离公式为：

，显然

，距离越靠近0，表示模式越相似；越靠近2，表明模式越不相似。把所有的模式距离加起来即时间序列的模式距离：

嗯……目前为止看起来还不错。但是需要注意的是，每一个模式可能跨越了不同的时间长度，而一个模式持续时间越长，它包含整个序列的信息就越多，这启发我们需要进行加权。因此：

，其中

，

为第i个模式所跨越的时间长度，

为总时间长度。

3. 形状距离shape distance

在模式距离的基础上，XIAO-LI DONG, CHENG-KUI GU, ZHENG-OU WANG提出了形状距离，进一步提升了度量效果。该算法于2006年在国际机器学习与控制会议上提出。

如果理解了pattern distance的计算，那么理解shape distance将会非常简单，因为shape distance归根到底、总而言之就是加了一个振幅的改变量并重新设定了模式序列。

假设我们已经得到了等模式数化后的序列：

设振幅amplitude改变量序列为A，则：

，就是我们每一个分割区间的端点对应的序列值值差，那么我们可以得到：

形状距离的最终计算公式为：

，同样

为时间权重。注意原序列S有N个点，模式序列和振幅改变量序列都是只有N-1个点。而这里的模式m也要重新定义。

下面以股票为例进行说明。我们认为股票的价格走势通常有七种状态：{加速下降，水平下降，减速下降，不变，减速上升，水平上升，加速上升}，我们用模式

来描述这一点。

现在我们设定一个阈值th来帮助我们区分这7种状态，设Ki表示通过PLR划分后的第K段直线的斜率，还记得在模式距离中怎么设定模式的吗？

如果斜率小于0，则表示下降，记为-1；斜率大于0，则表示上升，记为1；如果斜率等于0，则表示不变，那么记为0。

现在情况稍微复杂了一些，因为我们引入了更多的模式(7种)。 此时属于{-3，-2，-1}中的一种，那如何来知晓这支股票是加速下降、水平下降还是减速下降呢？斜率的改变量可以帮助我们，如果下一期斜率的改变量小于0，那说明斜率在递减，直线变得更陡峭了，是加速下降的，因此设为-3模式。如果

，那么是水平下降的，设为-2。如果

,说明是减速下降的，设为-1。

如果斜率的在这个范围内，那么就认为是近似不变的，设为0。

以此类推。

用一张表来总结一下：

形状距离公式是满足以下四个距离定理的：

此外，为了提高模型的准确度，较少绝对值的差异对整个模型的影响，一般需要对原序列进行标准化处理。完整的算法流程图如下：

形状距离由于加入了更多的模式，使得距离的度量更加精确，效果会好于模式距离。

4. DTW

讲了这么多距离，有没有更简单的方法来度量时间序列的相似性呢？相关系数可以吗？

相关系数是一个非常优美的指标，不仅能衡量相关性，也能衡量相似性，但是在时间序列中一般不使用相关系数衡量相似性，这是因为它不能解决图形平移的问题。

还是举个栗子~

y1与y2是平移的关系，显然与y1最相似的应该是y2，然而相关系数告诉我们是y3。

下面是Python进行模拟的代码：

import numpy as np

x=np.arange(-np.pi*2,np.pi*2,0.1)

y1=np.sin(x)

y2=np.cos(x)

y3=x

plt.plot(y1)

plt.plot(y2)

plt.plot(y3)

plt.legend(['y1','y2','y3'])

np.corrcoef([y1,y2,y3])

Out[57]:

array([[ 1.00000000e+00, -1.76804105e-04, -3.88194383e-01],

[-1.76804105e-04, 1.00000000e+00, -1.28528471e-02],

[-3.88194383e-01, -1.28528471e-02, 1.00000000e+00]])

DTW(dynamic time warping)动态时间规整就是专为解决此问题而生~

DTW实际上是计算欧式距离，我们在之前讲过，欧式距离不能很好地度量形状相似性。左边这幅图更清楚地表现了这一点，按照一般欧式距离的计算方法，所有点直接对应下来。而DTW就是找到序列之间正确对应的点再计算他们的距离。而由于DTW这一出色的特质，这使得DTW在语音识别领域有着广泛的应用。因为一个音节，它可能拖的很长或很短，如何去正确地识别相似声音或音节对于语音识别至关重要。图源自：https://blog.csdn.net/zouxy09/article/details/9140207

如上图所示，a对应的点应该是b而不是b'，DTW的根本任务就是将点进行正确的对应。

那正确对应的标准是什么呢？

DTW认为，如果两个序列的点正确对应了，那么他们的距离(欧式距离)达到最小。

OK，现在问题很简单了，就找距离最小的点。一个序列的一个点可能对应另一个序列的多个点，如果穷举出所有的可能去找出最合适的点，这在时间复杂度上属于一个NP难的问题。我们需要一个行之有效的算法——动态规划。

由于DTW的理论比较晦涩，直接看公式让人云里雾里，为了更清楚简单地把DTW讲明白，我用Python进行一遍模拟的计算。

简单点，说话的方式简单点~~~

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sns

x = np.array([1, 1, 2, 3, 2, 0])

y = np.array([0, 1, 1, 2, 3, 2, 1])

plt.plot(x,'r', label='x')

plt.plot(y, 'g', label='y')

plt.legend()

生成了两个长度不等的序列x,y。

distances = np.zeros((len(y), len(x)))

for i in range(len(y)):

for j in range(len(x)):

distances[i,j] = (x[j]-y[i])**2

distances

Out[60]:

array([[1., 1., 4., 9., 4., 0.],

[0., 0., 1., 4., 1., 1.],

[0., 0., 1., 4., 1., 1.],

[1., 1., 0., 1., 0., 4.],

[4., 4., 1., 0., 1., 9.],

[1., 1., 0., 1., 0., 4.],

[0., 0., 1., 4., 1., 1.]])

计算两个序列的距离矩阵。横着表示x序列，竖着是y序列，比如说第0行第0个元素1表示x序列的第0个值和y序列的第0个值的距离(Python的索引从0开始)，即

;再比如第4行第1列的元素4(如果你的第4行第1列元素没有找到4的话，请把索引从0开始)表示x序列第1个值和y序列的第4个值间的距离，即

。

我们来做个可视化，颜色越深表示距离越远，比如说最远的距离是x的第三个值和y的第0个值之间的距离。

这里注意一下，我们的图形和距离矩阵是没有一一对应的。我们的距离矩阵中，索引0从左上角开始，图形从左下角开始。

def distance_cost_plot(distances):

plt.imshow(distances, interpolation='nearest', cmap='Reds')

plt.gca().invert_yaxis()#倒转y轴，让它与x轴的都从左下角开始

plt.xlabel("X")

plt.ylabel("Y")

# plt.grid()

plt.colorbar()

distance_cost_plot(distances)

现在我们计算一个累积距离矩阵，累积距离最小是我们的最终目标。

accumulated_cost = np.zeros((len(y), len(x)))

accumulated_cost[0,0] = distances[0,0]

accumulated_cost

Out[80]:

array([[1., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.],

[0., 0., 0., 0., 0., 0.]])

distance_cost_plot(accumulated_cost)

显然累积距离矩阵的第0行第0列=距离矩阵的第0行第0列=1，我们必须经过起点吧……如果我们一直往右走，那么累积距离距离矩阵：

for i in range(1, len(x)):

accumulated_cost[0,i] = distances[0,i] + accumulated_cost[0, i-1]

accumulated_cost

array([[ 1., 2., 6., 15., 19., 19.],

[ 0., 0., 0., 0., 0., 0.],

[ 0., 0., 0., 0., 0., 0.],

[ 0., 0., 0., 0., 0., 0.],

[ 0., 0., 0., 0., 0., 0.],

[ 0., 0., 0., 0., 0., 0.],

[ 0., 0., 0., 0., 0., 0.]])

distance_cost_plot(accumulated_cost)

如果我们一直往上走，那么：

for i in range(1, len(y)):

accumulated_cost[i,0] = distances[i, 0] + accumulated_cost[i-1, 0]

accumulated_cost

array([[ 1., 2., 6., 15., 19., 19.],

[ 1., 0., 0., 0., 0., 0.],

[ 1., 0., 0., 0., 0., 0.],

[ 2., 0., 0., 0., 0., 0.],

[ 6., 0., 0., 0., 0., 0.],

[ 7., 0., 0., 0., 0., 0.],

[ 7., 0., 0., 0., 0., 0.]])

distance_cost_plot(accumulated_cost)

好了，累积距离矩阵的计算方式已经完全明白了。在继续下去之前，必须讲明累积距离矩阵中路径的意义。我们的目标是找一条路径使得累积距离最小。从(0,0)开始，经过这个方块表示将x,y的第0个点对应起来，当然，起点是必须经过的~

如果路径经过(1,0)，表示将x的第1个点和y序列的第0个点对应起来。

如果路径经过(1,1)，表示将x的第1个点和y序列的第1个点对应起来。

所以，路径的前进方向只有三个：向右，向上，斜右上。

这是显然，我们的路径不能后退，如果路径向右或向上了多个方块，这表明一个序列的一个点对应了另一个序列的多个点。

目前为止，我们已经完成了累积距离矩阵的两列的计算。对于任何其他的点，累积距离的新增只可能来自于三个方向：左边，下边，斜左下。所谓动态规划就是我每前进一步都选取使我当前行进距离最小的方向。因此，对于任何其他的点，累积距离的计算公式为：

现在，我们把累积距离矩阵计算完整：

for i in range(1, len(y)):

for j in range(1, len(x)):

accumulated_cost[i, j] = min(accumulated_cost[i-1, j-1], accumulated_cost[i-1, j], accumulated_cost[i, j-1]) + distances[i, j]

accumulated_cost

Out[86]:

array([[ 1., 2., 6., 15., 19., 19.],

[ 1., 1., 2., 6., 7., 8.],

[ 1., 1., 2., 6., 7., 8.],

[ 2., 2., 1., 2., 2., 6.],

[ 6., 6., 2., 1., 2., 11.],

[ 7., 7., 2., 2., 1., 5.],

[ 7., 7., 3., 6., 2., 2.]])

distance_cost_plot(accumulated_cost)

现在，最佳路径已经清晰地显示在了累积距离矩阵之中，就是图中颜色最淡的方块。

现在，我们只需要通过回溯的方法找回最佳路径就可以了：

path = [[len(x)-1, len(y)-1]]

i = len(y)-1

j = len(x)-1

while i>0 and j>0:

if i==0:

j = j - 1

elif j==0:

i = i - 1

else:

if accumulated_cost[i-1, j] == min(accumulated_cost[i-1, j-1], accumulated_cost[i-1, j], accumulated_cost[i, j-1]):

i = i - 1#来自于左边

elif accumulated_cost[i, j-1] == min(accumulated_cost[i-1, j-1], accumulated_cost[i-1, j], accumulated_cost[i, j-1]):

j = j-1#来自于下边

else:

i = i - 1#来自于左下边

j= j- 1

path.append([j, i])

path.append([0,0])

path

Out[89]: [[5, 6], [4, 5], [3, 4], [2, 3], [1, 2], [1, 1], [0, 1], [0, 0]]

path_x = [point[0] for point in path]

path_y = [point[1] for point in path]

distance_cost_plot(accumulated_cost)

plt.plot(path_x, path_y)

在刚才，我们已经讲明，最优路径所经过的地方表示两个序列应该对应的点。

计算这些点的欧式距离作为相似性的度量，这就是DTW。

Reference:Dong X L, Gu C K, Wang Z O. Research on shape-based time series similarity measure[C]//2006 International Conference on Machine Learning and Cybernetics. IEEE, 2006: 1253-1258.

Wang D, Rong G. Pattern distance of time series[J]. Zhejiang Daxue Xuebao(Gongxue Ban)/Journal of Zhejiang University(Engineering Science), 2004, 38(7): 795-798.

Kim S, Lee H, Ko H, et al. Pattern Matching Trading System Based on the Dynamic Time Warping Algorithm[J]. Sustainability, 2018, 10(12): 4641.

Keogh E, Chakrabarti K, Pazzani M, et al. Dimensionality reduction for fast similarity search in large time series databases[J]. Knowledge and information Systems, 2001, 3(3): 263-286.

Mori U, Mendiburu A, Lozano J A. Distance measures for time series in R: The TSdist package[J]. R journal, 2016, 8(2): 451-459.