时间序列相似性属于曲线相似性/曲线匹配(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在语音识别领域有着广泛的应用。因为一个音节，它可能拖的很长或很短，如何去正确地识别相似声音或音节对于语音识别至关重要。
如上图所示，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:
https://nipunbatra.github.io/blog/2014/dtw.html
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.
https://blog.csdn.net/xsdxs/article/details/86648605
知乎：如何判断两条轨迹（或曲线）的相似度？
张戎的博客
https://blog.csdn.net/qq_40006058/article/details/79992255
https://www.cnblogs.com/tornadomeet/archive/2012/03/23/2413363.html

时间序列相似性属于曲线相似性/曲线匹配(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.

面临的问题

当数据在时间线上不对齐的时候，使用传统的匹配方法，是无法使用传统的全局匹配度量法的。DTW是一种衡量两个时间序列之间的相似度的方法，主要应用在语音识别领域来识别两段语音是否表示同一个单词。

 

DTW原理

(Dynamic Time Warping, DTW) 动态时间归整

DTW通过把时间序列进行延伸和缩短，来计算两个时间序列性之间的相似性。

如下图所示，上下两条实线代表两个时间序列，时间序列之间的虚线代表两个时间序列之间的相似度，虚线的两端代表两个相似的点。



DTW使用所有这些相似点之间的距离的和，称之为归整路径距离(Warp Path Distance)来衡量两个时间序列之间的相似度。

给出两个序列：



Warping通常采用动态规划算法。为了对齐这两个序列，我们需要构造一个n x m的矩阵网格，矩阵元素(i, j)表示qi和cj两个点的距离d(qi, cj)，一般采用欧式距离，（也可以理解为失真度）。每一个矩阵元素(i, j)表示点qi和cj的对齐。

DP(dynamic programming)算法可以归结为寻找一条通过此网格中若干格点的路径，路径通过的格点即为两个序列进行计算的对齐的点。

有三个性质：

边界条件：w1=(1, 1)和wK=(m, n)。两个人分别说了同一个单词，但是由于语速、语气、语调等等各不相同，会导致采样得到的数据无法对齐。但是两段语音采样的第一个采样值和最后一个采样值肯定是两两对应的。
	
连续性：如果wk-1= (a', b')，那么对于路径的下一个点wk=(a, b)需要满足 (a-a') <=1和 (b-b') <=1。也就是不可能跨过某个点去匹配，只能和自己相邻的点对齐。这样可以保证Q和C中的每个坐标都在W中出现。
	
单调性：如果wk-1= (a', b')，那么对于路径的下一个点wk=(a, b)需要满足0<=(a-a’)和0<= (b-b’)。这限制W上面的点必须是随着时间单调进行的。以保证上图中的虚线不会相交。
由连续性和单调性可知，每次格点(i, j)前进方向只有三种：(i+1, j)，(i, j+1) 或 (i+1, j+1)。我们的目的是使得下面的规整代价最小的路径：



分母中的K主要是用来对不同的长度的规整路径做补偿。

这里我们定义一个累加距离(cumulative distances)。从(0, 0)点开始匹配这两个序列Q和C，每到一个点，之前所有的点计算的距离都会累加。到达终点(n, m)后，这个累积距离就是我们上面说的最后的总的距离，也就是序列Q和C的相似度。

示例：

对于两个序列：

X：3,5,6,7,7,1

Y：3,6,6,7,8,1,1

X和Y的距离矩阵如下


然后根据距离矩阵生成损失矩阵(Cost Matrix)或者叫累积距离矩阵 ，其计算方法如下：


	
第一行第一列元素为 的第一行第一列元素，在这里就是0；
	
	

	
其他位置的元素 的值则需要逐步计算，具体值的计算方法为 ，得到的如下：
	


 最后，两个序列的距离，由损失矩阵最后一行最后一列给出，在这里也就是2。

给定了距离矩阵，如何找到一条从左上角到右下角的路径，使得路径经过的元素值之和最小。

 

python 代码实现

 代码里用到的sys.maxsize，需要python3支持

import sys

def distance(a,b):
    return abs(a-b)

def Min(a,b,c):
    return min(a,min(b,c))

def dtw(X, Y):
    M = [[distance(X[j], Y[i]) for i in range(len(Y))] for j in range(len(X))]
    l1 = len(X)
    l2 = len(Y)
    D = [[0 for i in range(l2+1)] for i in range(l1+1)]
    D[0][0] = 0
    for i in range(1, l2 + 1):
        D[0][i] = sys.maxsize
    for j in range(1, l1 + 1):
        D[j][0] = sys.maxsize
    for j in range(1, l2+1):
        for i in range(1, l1+1):
            D[i][j] = M[i - 1][j - 1] + Min(D[i - 1][j], D[i][j - 1], D[i - 1][j - 1])
    return D[l1][l2]

if __name__ == '__main__':
    X = [3,5,6,7,7,1]
    Y = [3,6,6,7,8,1,1]
    Z = [2,5,7,7,7,7,2]
    value = dtw(X, Y)
    print(value)
    value = dtw(X, Z)
    print(value)


DTW虽然使用线性规划可以快速的求解，但是在面对比较长的时间序列是，O(N2)的时间复杂度还是很大。已经有很多改进的快速DTW算法，比如FastDTW等。

DTW常用加速手段

  常用的DTW加速手段:

(1). 限制。亦即减少累积距离矩阵的搜索空间，下图中阴影部分为实际的探索空间，空白的部分不进行探索。



 (2). 数据抽象。亦即把之前长度为N的时间序列规约成长度为M(M<N)表述方式：

 

 (3). 索引。索引是在进行分类和聚类时减少需要运行的DTW的次数的方法，并不能加速一次的DTW计算。

 FastDTW

  FastDTW综合使用限制和数据抽象两种方法来加速DTW的计算，主要分为三个步骤：

  (1). 粗粒度化。亦即首先对原始的时间序列进行数据抽象，数据抽象可以迭代执行多次1/1->1/2->1/4->1/16，粗粒度数据点是其对应的多个细粒度数据点的平均值。

  (2). 投影。在较粗粒度上对时间序列运行DTW算法。

  (3). 细粒度化。将在较粗粒度上得到的归整路径经过的方格进一步细粒度化到较细粒度的时间序列上。除了进行细粒度化之外，我们还额外的在较细粒度的空间内额外向外(横向，竖向，斜向)扩展K个粒度，K为半径参数，一般取为1或者2.

  FastDTW算法的具体执行流程如下图所示:



  第一个图表示在较粗粒度空间(1/8)内执行DTW算法。第二个图表示将较粗粒度空间(1/8)内求得的归整路径经过的方格细粒度化，并且向外(横向，竖向，斜向)扩展一个(由半径参数确定)细粒度单位后，再执行DTW得到的归整路径。第三个图和第四个图也是这样。

  由于采取了减少搜索空间的策略，FastDTW并不一定能够求得准确的DTW距离，但是FastDTW算法的时间复杂度比较低，为O(N)。

 

参考：

https://blog.csdn.net/raym0ndkwan/article/details/45614813

https://www.cnblogs.com/kemaswill/archive/2013/04/18/3028610.html

https://www.cnblogs.com/kemaswill/archive/2013/04/18/3029078.html