时间序列 相似 python 时间序列相似度检测

时间序列分析 - 23 DTW (时序相似度度量算法) 上

DTW初探

简介

    在时序分析中，DTW(Dynamic Time Warping)是用来检测两个时序相似程度的算法，而这个相似程度通常用一个距离来表示。例如如下的两个序列，

$a=[1,2,3]$

$b=[3,2,2]$

我们该如何衡量这两个序列的距离呢？一个比较明显的方法是对 𝑎 , 𝑏 这两个序列中的元素按照位置一一计算距离，最后加总或者加总后平均。听上去是一个比较简单其易于理解的方法。我们可以观察一下下图理解一下时序相似程度度量。

但如果是如下时序，

$a=[1,2,3]$

$b=[2,2,2,3,4]$

我们又该怎么办呢？

DTW就是一种最常用的解决此类问题的算法。

    前面所提到的简单算法实际上就是计算两个对齐向量的欧几里得距离或者曼哈顿距离，这种方法不但无法解决序列长度不一样问题，而且即使序列长度一样时得到的距离也不是一个好的度量。

简单来说，DTW算法提供了两个序列的具有弹性的对齐方式，寻找最佳对齐方式，并基于此最佳对齐方式再计算距离。

应用

    DTW有非常广泛的应用。例如它可以用来发现相似的步行方式，即使一个走得快另一个走得慢，亦或一个在减速或加速。

DTW已经被应用到了视频、音频、图像等多种可以序列化的数据，例如根据音频识别说话的人、签名识别等。

    另外一个典型应用是在股市的时序分析上。在股票市场中，人们总是希望可以预测未来，但是一般的机器学习算法在这个领域并不是非常有效。因为机器学习的预测任务一般需要其训练集和测试集拥有同样维度的特征属性。但是我们知道即便是一个股票的非常相同的模式出现了，它们的K线和其他指标的反映长度也是不大相同的。

DTW 算法规则

DTW是在一定的规则和约束条件下，计算两个序列的最优匹配：

• 第一个序列的每个位置元素必须匹配第二个序列的一个或多个位置元素，反之亦然。
• 第一个序列的第一个元素必须匹配第二个序列的第一个元素，但可以还匹配其他位置。
• 第一个序列的最后一个元素必须匹配第二个序列的最后一个元素，但可以还匹配其他位置。
• 两个序列的位置映射关系必须是单调增的，反之亦然。举例来说就是假如第一个序列的1，3位置映射第二个序列的2，3位置是可以的，如果映射到3，2位置就是不行的。

最优匹配就是满足上面的所有约束规则之下，具备最小损失的匹配。其中最小损失就是映射位置上的数值的距离的加总。 总结一下就是头尾必须匹配，没有交叉映射，没有残余子序列。

先来看一下伪码：  
int DTWDistance( s: array [1..n],  t: array [1..m] ) {  
    DTW := array [0..n, 0..m]
    
    for i := 0 to n
        for j := 0 to m
            DTW[i, j] := infinity
    DTW[0, 0] := 0
    
    for i := 1 to n
        for j := 1 to m
            cost := d(s[i], t[j])
            DTW[i, j] := cost + minimum(DTW[i-1, j  ],    // insertion
                                        DTW[i  , j-1],    // deletion
                                        DTW[i-1, j-1])    // match
    
    return DTW[n, m]
}

这里，$d(x,y)=|x-y|$，而 $DTW[i,j]$ 代表在最好的对齐之下，序列 $s[1:i]$ 与 $t[1:j]$ 之间的距离。
这个算法的核心是:

DTW[i, j] := cost + minimum(DTW[i-1, j ], //insertion
                   DTW[i  , j-1],  //deletion
                   DTW[i-1, j-1])  //match

很明显这是一个动态规划算法。两个序列$s[1:i]$与$t[1:j]$之间的距离就是$i$与$j$的距离加上下面三个子问题的最小距离

• 子问题1：DTW[i-1,j] 理解为$s$序列插入一个点
• 子问题2：DTW[i,j-1] 理解为$t$序列删除一个点
• 子问题3：DTW[i-1,j-1] 理解为两个序列在上一个时间点匹配

Python代码如下

def dtw(s, t):
    n, m = len(s), len(t)
    dtw_matrix = np.zeros((n+1, m+1))
    for i in range(n+1):
        for j in range(m+1):
            dtw_matrix[i, j] = np.inf
    dtw_matrix[0, 0] = 0
    
    for i in range(1, n+1):
        for j in range(1, m+1):
            cost = abs(s[i-1] - t[j-1])
            # take last min from a square box
            last_min = np.min([dtw_matrix[i-1, j], dtw_matrix[i, j-1], dtw_matrix[i-1, j-1]])
            dtw_matrix[i, j] = cost + last_min
    return dtw_matrix

举个例子

import numpy as np
a = [1,2,3]
b = [2,2,2,3,4]
dtw(a,b)

array([[ 0., inf, inf, inf, inf, inf],
[inf, 1., 2., 3., 5., 8.],
[inf, 1., 1., 1., 2., 4.],
[inf, 2., 2., 2., 1., 2.]])
最后一个元素的值就是两个序列的距离，为2。

增加窗口限制

上面的算法存在一个问题，就是允许一个序列中的一个元素匹配另一个序列中的元素个数是没有限制的。在这种情况下，如果其中一个序列很长的化将会导致性能问题。例如，
$a=[1,2,3]$
$b=[1,2,2,2,2,2,...,5]$
为了避免这个问题，我们可以增加一个窗口参数来限制一个时序点可以在另一个时序匹配的元素的个数。

伪码如下

int DTWDistance( s: array [1..n], t: array [1..m], w: int ) {   
    
    DTW := array [0..n, 0..m]

    w := max(w, abs(n-m)) // adapt window size (*)

    for i := 0 to n
        for j:= 0 to m
            DTW[i, j] := infinity
    DTW[0, 0] := 0
    for i := 1 to n
        for j := max(1, i-w) to min(m, i+w)
            DTW[i, j] := 0

    for i := 1 to n
        for j := max(1, i-w) to min(m, i+w)
            cost := d(s[i], t[j])
            DTW[i, j] := cost + minimum(DTW[i-1, j  ],    // insertion
                               DTW[i  , j-1],    // deletion
                               DTW[i-1, j-1])    // match

    return DTW[n, m]
}

Python代码如下

def dtw(s, t, window):
    n, m = len(s), len(t)
    w = np.max([window, abs(n-m)])
    dtw_matrix = np.zeros((n+1, m+1))
    
    for i in range(n+1):
        for j in range(m+1):
            dtw_matrix[i, j] = np.inf
    dtw_matrix[0, 0] = 0
    
    for i in range(1, n+1):
        for j in range(np.max([1, i-w]), np.min([m, i+w])+1):
            dtw_matrix[i, j] = 0
    
    for i in range(1, n+1):
        for j in range(np.max([1, i-w]), np.min([m, i+w])+1):
            cost = abs(s[i-1] - t[j-1])
            # take last min from a square box
            last_min = np.min([dtw_matrix[i-1, j], dtw_matrix[i, j-1], dtw_matrix[i-1, j-1]])
            dtw_matrix[i, j] = cost + last_min
    return dtw_matrix

a = [1,2,3,3,5]
b = [1,2,2,2,2,2,2,4]
dtw(a,b,window=3)

array([[ 0., inf, inf, inf, inf, inf, inf, inf, inf],
 [inf, 0., 1., 2., 3., inf, inf, inf, inf],
 [inf, 1., 0., 0., 0., 0., inf, inf, inf],
 [inf, 3., 1., 1., 1., 1., 1., inf, inf],
 [inf, 5., 2., 2., 2., 2., 2., 2., inf],
 [inf, inf, 5., 5., 5., 5., 5., 5., 3.]])

在带有窗口参数控制下，结果为3.

Python 包

前面的示例代码采用了递归来实现动态规划，这种方式虽然实现简单但是运行效率很低。所以在实际的场景中一般会采用一些优化算法，例如
PrunedDTW, SparseDTW, FastDTW 和 MultiscaleDTW 这里我们采用FastDTW来试演一下

from fastdtw import fastdtw
from scipy.spatial.distance import euclidean

x = np.array([1, 2, 3, 3, 7])
y = np.array([1, 2, 2, 2, 2, 2, 2, 4])

distance, path = fastdtw(x, y, dist=euclidean)

print(distance)
print(path)

# 5.0
# [(0, 0), (1, 1), (1, 2), (1, 3), (1, 4), (2, 5), (3, 6), (4, 7)]

5.0
[(0, 0), (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (2, 7), (3, 7), (4, 7)]

未完，待续…