【ML-7-應用】聚類算法-時間序列聚類（DTW和LB_Keogh距離）

目錄

1. 問題分析
2. 數據處理
3. 代碼實現
4. 結果

今天兄弟找我幫忙：希望基於白細胞數據把新型肺炎患者的數據做一個聚類並畫出大體曲線：也就是將相同變化的曲線進行分類並擬合。定位此問題為無監督的分類問題。因此想到了聚類的方法。

一、問題分析

1、首先嘗試了使用：提取時間序列的統計學特征值，例如最大值，最小值等。然后利目前常用的算法根據提取的特征進行分類，例如Naive Bayes, SVMs，KNN 等。發現效果並不是很好。

2、嘗試基於K-means的無監督形式分類，這種分類方式基於兩個數據的距離進行分類，這樣要定義號距離的概念，后來查閱資料，考慮使用動態時間規整（Dynamic Time Warping, DTW）。下文主要基於這方面進行展開。

二、數據處理

給出的數據較為完整，就一個excel表格，做了以下簡單的排序，原始數據可見文末github地址。

三、代碼實現

3.1 動態時間規整（Dynamic Time Warping, DTW）

如果是歐拉距離：則ts3比ts2更接近ts1，但是肉眼看並非如此。故引出DTW距離。

  動態時間規整算法，故名思議，就是把兩個代表同一個類型的事物的不同長度序列進行時間上的"對齊"。比如DTW最常用的地方，語音識別中，同一個字母，由不同人發音，長短肯定不一樣，把聲音記錄下來以后，它的信號肯定是很相似的，只是在時間上不太對整齊而已。所以我們需要用一個函數拉長或者縮短其中一個信號，使得它們之間的誤差達到最小。下面這篇博文給了比較好的解釋：https://blog.csdn.net/lin_limin/article/details/81241058。 簡單英文解釋如下（簡而言之：就是允許錯開求差值，並且取最小的那個作為距離。）

DTW距離代碼定義如下：

    1 	def DTWDistance(s1, s2):
    2 	    DTW={}
    3 	    for i in range(len(s1)):
    4 	        DTW[(i, -1)] = float('inf')
    5 	    for i in range(len(s2)):
    6 	        DTW[(-1, i)] = float('inf')
    7 	    DTW[(-1, -1)] = 0
    8 	    for i in range(len(s1)):
    9 	        for j in range(len(s2)):
   10 	            dist= (s1[i]-s2[j])**2
   11 	            DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], DTW[(i-1, j-1)])
   12 			
   13 	    return math.sqrt(DTW[len(s1)-1, len(s2)-1])
   14 
   15 

這樣求解相對較為麻煩，時間復雜度比較高，故做了一個小的加速：

#DTW距離,只檢測前W個窗口的值，定義錯開的部分W，減少遞歸尋找量

    1 	def DTWDistance(s1, s2, w):
    2 	    DTW = {}
    3 	    w = max(w, abs(len(s1) - len(s2)))
    4 	    for i in range(-1, len(s1)):
    5 	        for j in range(-1, len(s2)):
    6 	            DTW[(i, j)] = float('inf')
    7 	    DTW[(-1, -1)] = 0
    8 	    for i in range(len(s1)):
    9 	        for j in range(max(0, i - w), min(len(s2), i + w)):
   10 	            dist = (s1[i] - s2[j]) ** 2
   11 	            DTW[(i, j)] = dist + min(DTW[(i - 1, j)], DTW[(i, j - 1)], DTW[(i - 1, j - 1)])
   12 	    return math.sqrt(DTW[len(s1) - 1, len(s2) - 1])

3.2 LB_Keogh距離

主要思想是在搜索數據很大的時候， 逐個用DTW算法比較每一條是否匹配非常耗時。那我們能不能使用一種計算較快的近似方法計算LB， 通過LB處理掉大部分不可能是最優匹配序列的序列，對於剩下的序列在使用DTW逐個比較呢？英文解釋如下：

中文解釋如下，主要參考其它博文：LB_keogh的定義相對復雜，包括兩部分。

第一部分為Q的{U， L} 包絡曲線（具體如圖）， 給Q序列的每個時間步定義上下界。 定義如下

其中 r 是一段滑行窗距離，可以自定義。

示意圖如下：

U 為上包絡線，就是把每個時間步為Q當前時間步前后r的范圍內最大的數。L 下包絡線同理。那么LB_Keogh定義如下：

用圖像描述如下：

    1 	def LB_Keogh(s1, s2, r):
    2 	    LB_sum = 0
    3 	    for ind, i in enumerate(s1):
    4 	        # print(s2)
    5 	        lower_bound = min(s2[(ind - r if ind - r >= 0 else 0):(ind + r)])
    6 	        upper_bound = max(s2[(ind - r if ind - r >= 0 else 0):(ind + r)])
    7 	        if i >= upper_bound:
    8 	            LB_sum = LB_sum + (i - upper_bound) ** 2
    9 	        elif i < lower_bound:
   10 	            LB_sum = LB_sum + (i - lower_bound) ** 2
   11 	    return math.sqrt(LB_sum)

3.3 使用k-means算法實現聚類

    1 	#3  定義K-means算法
    2 	#num_clust分類的數量，
    3 	def k_means_clust(data, num_clust, num_iter, w=3):
    4 	    ## 步驟一: 初始化均值點
    5 	    centroids = random.sample(list(data), num_clust)
    6 	    counter = 0
    7 	    for n in range(num_iter):
    8 	        counter += 1
    9 	        # print
   10 	        # counter
   11 	        assignments = {}  #存儲類別0，1，2等類號和所包含的類的號碼
   12 	        # 遍歷每一個樣本點 i ,因為本題與之前有所不同，多了ind的編碼
   13 	        for ind, i in enumerate(data):
   14 	            min_dist = float('inf')   #最近距離，初始定一個較大的值
   15 	            closest_clust = None     # closest_clust：最近的均值點編號
   16 	            ## 步驟二: 尋找最近的均值點
   17 	            for c_ind, j in enumerate(centroids):  #每個點和中心點的距離，共有num_clust個值
   18 	                if LB_Keogh(i, j, 3) < min_dist:    #循環去找最小的那個
   19 	                    cur_dist = DTWDistance(i, j, w)
   20 	                    if cur_dist < min_dist:         #找到了ind點距離c_ind最近
   21 	                        min_dist = cur_dist
   22 	                        closest_clust = c_ind
   23 	            ## 步驟三: 更新 ind 所屬簇
   24 	            # print(closest_clust)
   25 	            if closest_clust in assignments:
   26 	                assignments[closest_clust].append(ind)
   27 	            else:
   28 	                assignments[closest_clust] = []
   29 	                assignments[closest_clust].append(ind)
   30 	        # recalculate centroids of clusters  ## 步驟四: 更新簇的均值點
   31 	        for key in assignments:
   32 	            clust_sum = 0
   33 	            for k in assignments[key]:
   34 	                clust_sum = clust_sum + data[k]
   35 	            centroids[key] = [m / len(assignments[key]) for m in clust_sum]
   36 	    return centroids,assignments    #返回聚類中心值，和聚類的所有點的數組序號

3.4 根據聚類打印出具體分類情況：

    1 	num_clust = 2  #定義需要分類的數量
    2 	centroids,assignments = k_means_clust(WBCData, num_clust,800, 3)
    3 	for i in range(num_clust):
    4 	    s = []
    5 	    WBC01 = []
    6 	    days01 = []
    7 	    for j, indj in enumerate(assignments[i]):  #畫出各分類點的坐標
    8 	        s.append(int(Numb[indj*30]))
    9 	        WBC01 = np.hstack((WBC01,WBC[30 * indj:30 * indj + 30]))
   10 	        days01 = np.hstack((days01 , days[0: 30]))
   11 	    print(s)
   12 	    plt.title('%s' % s)
   13 	    plt.plot(centroids[i],c="r",lw=4)
   14 	    plt.scatter(days01, WBC01 )
   15 	    plt.show()

四、結果

定義了分成兩類的情形，可以根據num_clust 的值進行靈活的調整，等於2是的分類和圖示情況如下：

WBC01：[6774, 7193, 8070, 8108, 8195, 2020006799, 2020007003, 2020007251, 2020007420, 2020007636, 2020007718, 2020007928, 2020007934, 2020008022, 2020008196, 2020008239, ……] 不全部列出

WBC02：[2020007250, 2020007388, 2020007389, 2020007422, 2020007625, 2020007703, 2020007927 , ……]

說明：

代碼訓練過程中，一定要注意數據類型，比如matrix和ndarray,雖然打印的時候都是（45，30），但是再訓練的時候，稍加不注意，就會導致亂七八糟的問題，需要打印排查好久。

本文的數據和代碼，請登錄：my github ,進行下載。若是對您有用，請不吝給顆星。

附件一：整體代碼

    1 	# Author:yifan
    2 	
    3 	import pandas as pd
    4 	import numpy as np
    5 	import matplotlib.pylab as plt
    6 	import math
    7 	import random
    8 	import  xlrd
    9 	import  numpy as np
   10 	import sys
   11 	
   12 	
   13 	#演示
   14 	# import pandas as pd
   15 	# import numpy as np
   16 	# import matplotlib.pylab as plt
   17 	# x=np.linspace(0,50,100)
   18 	# ts1=pd.Series(3.1*np.sin(x/1.5)+3.5)
   19 	# ts2=pd.Series(2.2*np.sin(x/3.5+2.4)+3.2)
   20 	# ts3=pd.Series(0.04*x+3.0)
   21 	# ts1.plot()
   22 	# ts2.plot()
   23 	# ts3.plot()
   24 	# plt.ylim(-2,10)
   25 	# plt.legend(['ts1','ts2','ts3'])
   26 	# plt.show()
   27 	# def euclid_dist(t1,t2):
   28 	#     return math.sqrt(sum((t1-t2)**2))
   29 	# print (euclid_dist(ts1,ts2))   #26.959216038
   30 	# print (euclid_dist(ts1,ts3))   #23.1892491903
   31 	
   32 	
   33 	#1   數據提取
   34 	# workbook = xlrd.open_workbook(r"D:\datatest\xinguanfeiyan\20200229.xlsx")
   35 	workbook = xlrd.open_workbook(r"D:\datatest\xinguanfeiyan\20200315L.xlsx")
   36 	worksheet = workbook.sheet_by_index(0)
   37 	n = worksheet.nrows
   38 	days = [0]*(n-1)
   39 	WBC = [0]*(n-1)
   40 	Numb = [0]*(n-1)
   41 	for i in range(1,n):
   42 	    Numb[i-1] = worksheet.cell_value(i,0)
   43 	    days[i-1] = worksheet.cell_value(i,1)
   44 	    WBC[i-1] = worksheet.cell_value(i,2)
   45 	s=int(len(WBC)/30)
   46 	WBCData = np.mat(WBC).reshape(s,30)
   47 	WBCData = np.array(WBCData)
   48 	# print(WBCData)
   49 	
   50 	#2  定義相似距離
   51 	#DTW距離，時間復雜度為兩個時間序列長度相乘
   52 	def DTWDistance(s1, s2):
   53 	    DTW = {}
   54 	    for i in range(len(s1)):
   55 	        DTW[(i, -1)] = float('inf')
   56 	    for i in range(len(s2)):
   57 	        DTW[(-1, i)] = float('inf')
   58 	    DTW[(-1, -1)] = 0
   59 	    for i in range(len(s1)):
   60 	        for j in range(len(s2)):
   61 	            dist = (s1[i] - s2[j]) ** 2
   62 	            DTW[(i, j)] = dist + min(DTW[(i - 1, j)], DTW[(i, j - 1)], DTW[(i - 1, j - 1)])
   63 	    return math.sqrt(DTW[len(s1) - 1, len(s2) - 1])
   64 	
   65 	#DTW距離,只檢測前W個窗口的值，減少時間復雜度
   66 	def DTWDistance(s1, s2, w):
   67 	    DTW = {}
   68 	    w = max(w, abs(len(s1) - len(s2)))
   69 	    for i in range(-1, len(s1)):
   70 	        for j in range(-1, len(s2)):
   71 	            DTW[(i, j)] = float('inf')
   72 	    DTW[(-1, -1)] = 0
   73 	    for i in range(len(s1)):
   74 	        for j in range(max(0, i - w), min(len(s2), i + w)):
   75 	            dist = (s1[i] - s2[j]) ** 2
   76 	            DTW[(i, j)] = dist + min(DTW[(i - 1, j)], DTW[(i, j - 1)], DTW[(i - 1, j - 1)])
   77 	    return math.sqrt(DTW[len(s1) - 1, len(s2) - 1])
   78 	#Another way to speed things up is to use the LB Keogh lower bound of dynamic time warping
   79 	def LB_Keogh(s1, s2, r):
   80 	    LB_sum = 0
   81 	    for ind, i in enumerate(s1):
   82 	        # print(s2)
   83 	        lower_bound = min(s2[(ind - r if ind - r >= 0 else 0):(ind + r)])
   84 	        upper_bound = max(s2[(ind - r if ind - r >= 0 else 0):(ind + r)])
   85 	        if i >= upper_bound:
   86 	            LB_sum = LB_sum + (i - upper_bound) ** 2
   87 	        elif i < lower_bound:
   88 	            LB_sum = LB_sum + (i - lower_bound) ** 2
   89 	    return math.sqrt(LB_sum)
   90 	
   91 	
   92 	#3  定義K-means算法
   93 	#num_clust分類的數量，
   94 	def k_means_clust(data, num_clust, num_iter, w=3):
   95 	    ## 步驟一: 初始化均值點
   96 	    centroids = random.sample(list(data), num_clust)
   97 	    counter = 0
   98 	    for n in range(num_iter):
   99 	        counter += 1
  100 	        # print
  101 	        # counter
  102 	        assignments = {}  #存儲類別0，1，2等類號和所包含的類的號碼
  103 	        # 遍歷每一個樣本點 i ,因為本題與之前有所不同，多了ind的編碼
  104 	        for ind, i in enumerate(data):
  105 	            min_dist = float('inf')   #最近距離，初始定一個較大的值
  106 	            closest_clust = None     # closest_clust：最近的均值點編號
  107 	            ## 步驟二: 尋找最近的均值點
  108 	            for c_ind, j in enumerate(centroids):  #每個點和中心點的距離，共有num_clust個值
  109 	                if LB_Keogh(i, j, 3) < min_dist:    #循環去找最小的那個
  110 	                    cur_dist = DTWDistance(i, j, w)
  111 	                    if cur_dist < min_dist:         #找到了ind點距離c_ind最近
  112 	                        min_dist = cur_dist
  113 	                        closest_clust = c_ind
  114 	            ## 步驟三: 更新 ind 所屬簇
  115 	            # print(closest_clust)
  116 	            if closest_clust in assignments:
  117 	                assignments[closest_clust].append(ind)
  118 	            else:
  119 	                assignments[closest_clust] = []
  120 	                assignments[closest_clust].append(ind)
  121 	        # recalculate centroids of clusters  ## 步驟四: 更新簇的均值點
  122 	        for key in assignments:
  123 	            clust_sum = 0
  124 	            for k in assignments[key]:
  125 	                clust_sum = clust_sum + data[k]
  126 	            centroids[key] = [m / len(assignments[key]) for m in clust_sum]
  127 	    return centroids,assignments
  128 	
  129 	#
  130 	num_clust = 2  #定義需要分類的數量
  131 	centroids,assignments = k_means_clust(WBCData, num_clust,800, 3)
  132 	for i in range(num_clust):
  133 	    s = []
  134 	    WBC01 = []
  135 	    days01 = []
  136 	    for j, indj in enumerate(assignments[i]):  #畫出各分類點的坐標
  137 	        s.append(int(Numb[indj*30]))
  138 	        WBC01 = np.hstack((WBC01,WBC[30 * indj:30 * indj + 30]))
  139 	        days01 = np.hstack((days01 , days[0: 30]))
  140 	    print(s)
  141 	    plt.title('%s' % s)
  142 	    plt.plot(centroids[i],c="r",lw=4)
  143 	    plt.scatter(days01, WBC01 )
  144 	    plt.show()