異常檢測(1)——局部異常因子算法

時間 2019-12-13

標籤異常檢測局部因子算法简体版

原文原文鏈接

　　局部異常因子算法（Local Outlier Factor）經過計算「局部可達密度」來反映一個樣本的異常程度，一個樣本點的局部可達密度越大，這個點就越有多是異常點。html

k距離和k距離鄰域

　　某一點P的k距離（k-distance）很容易解釋，就是點P和距離點P第k近的點之間距離，但不包括P。假設P是學校，葛小倫、劉闖、趙信、薔薇、琪琳、炙心6個同窗都住在學校附近：算法

　　圖1緩存

　　爲了簡單前起見，將P放置在原點。用歐幾里得距離表示每一個同窗的家到學校的距離，由近及遠分別是葛小倫<劉闖<趙信<薔薇=炙心<琪琳，距P最近的是葛小倫，第2近的是劉闖，第3近的是趙信，第4近的是薔薇和炙心，第5近的是琪琳。app

　　當k=t時，在數據集中用下式表示x⁽ⁱ⁾的k距離，其中x^(k=t)表示距x⁽ⁱ⁾第k遠的數據樣本：dom

　　所謂P的k距離鄰域（k-distance neighborhood of p），就是到P的直線距離小於P的k距離的全部數據樣本構成的集合。在圖9.1中，當k=3時，P的k距離鄰域是{葛小倫、劉闖、趙信}；當k=4時，P的k距離鄰域是{葛小倫、劉闖、薔薇、炙心}。post

　　能夠把P看做圓心，把P的k距離看做半徑作一個圓，圓中的樣本點就是P的k距離鄰域：學習

　　圖2測試

可達距離

　　x⁽ⁱ⁾到x^(j)的可達距離（Rechabiliby Distance）能夠表示爲：spa

　　它的含義是，當x⁽ⁱ⁾距x^(j)的距離比x⁽ⁱ⁾距x^(k=t)更近時，直接用數值較大的|| x⁽ⁱ⁾- x^(k=k)||表示x⁽ⁱ⁾到x^(j)的可達距離，不然用|| x⁽ⁱ⁾- x^(j)||表示。3d

　　以圖1爲例，當k=3時，距P第3近的同窗是趙信，P的k距離鄰域是{葛小倫、劉闖、趙信}，則P到他們的可達距離都等於學校到趙信家的距離；P到薔薇、炙心、琪琳的可達距離與學校到她們的實際距離相等：

　　因爲以x⁽ⁱ⁾的k近距離和x^(j)的k距離並不相等，因此x⁽ⁱ⁾到x^(j)的可達距離和x^(j)到x⁽ⁱ⁾的可達距離並不相等：

　　圖3

　　在圖3中，當k=2時，先計算x⁽¹⁾到x⁽²⁾的可達距離，此時距x⁽¹⁾第2近的點是x⁽⁴⁾：

　　再來看看x⁽²⁾到x⁽¹⁾的可達距離，距x⁽²⁾第2近的點是x⁽³⁾：

　　因而可知：

局部可達密度

　　點與點之間的密度很容易理解，點之間距離越遠，密度越低，距離越近，密度越高。局部可達密度與整體密度相似，只不過是用k距離鄰域計算的，因此稱爲「局部」。

　　假設x⁽ⁱ⁾的k距離鄰域中有N個樣本點，用x_N表示，x_N^(t)表示x_N中的第t個樣本點，那麼x⁽ⁱ⁾的局部可達密度（local reachability density）能夠寫成：

　　上式中x⁽ⁱ⁾的局部可達密度是x⁽ⁱ⁾的k距離鄰域中全部樣本點到x⁽ⁱ⁾的可達距離的平均值的倒數。表明x_N中樣本點的密集程度，密集程度越高，該值越小，它的倒數（x⁽ⁱ⁾的局部可達密度）的值越大，x⁽ⁱ⁾和x_N越多是同一簇；反之，若是x⁽ⁱ⁾是異常數據，那麼x_N中的樣本點到x⁽ⁱ⁾的可達距離將會取這些樣本點到x⁽ⁱ⁾的直線距離，這個距離將遠大於x_N中樣本點的k距離，最終致使較大，它的倒數（x⁽ⁱ⁾的局部可達密度）的值較小。簡單而言，x⁽ⁱ⁾的局部可達密度越大，x⁽ⁱ⁾越靠近它鄰域中的點；x⁽ⁱ⁾的局部可達密度越小，x⁽ⁱ⁾越遠離它鄰域中的點。

局部異常因子

　　有了局部可達密度，就能夠進一步定義局部異常因子（local outlier factor）：

　　上式的分子表示x⁽ⁱ⁾的k距離鄰域中的全部樣本點的局部可達密度的均值，分母是x⁽ⁱ⁾的局部可達密度。它其實是經過比較x⁽ⁱ⁾的密度和其鄰域的密度來判斷x⁽ⁱ⁾是不是異常點，x⁽ⁱ⁾的密度越低，LRD_k(x⁽ⁱ⁾)越小，LOF_k(x⁽ⁱ⁾)的值越大，x⁽ⁱ⁾越多是異常點；x⁽ⁱ⁾的密度越高，LRD_k(x⁽ⁱ⁾)越大，LOF_k(x⁽ⁱ⁾)的值越接近1或小於1，x⁽ⁱ⁾越多是正常的樣本點。

算法實現

　　咱們已經介紹了局部異常因子算法，來看看它是如何工做的。先來建立一些測試數據：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 
 4 def create_train():
 5     np.random.seed(42) # 設置seed使每次生成的隨機數都相等
 6     # 生成100個2維數據，它們是以0爲均值、以1爲標準差的正態分佈
 7     X_inliers = 0.3 * np.random.randn(100, 2)
 8     # 構造兩組間隔必定距離的樣本點做爲訓練數據
 9     X_inliers = np.r_[X_inliers + 2, X_inliers - 2]
10     # 構造20個可能的異常數據，從一個均勻分佈[low,high)中隨機採樣
11     X_outliers = np.random.uniform(low=-4, high=4, size=(20, 2))
12     # 將X_inliers和X_outliers鏈接起來做爲訓練集
13     return np.r_[X_inliers, X_outliers]
14 
15 def display(X):
16     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
17     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
18     plt.title("局部異常因子(LOF)")
19     plt.scatter(X[:, 0], X[:, 1], color='k', s=3., label='樣本點')
20 
21     plt.axis('tight') # 讓座標軸的比例尺適應數據量
22     plt.xlim((-5, 5)) # x座標軸的刻度範圍
23     plt.ylim((-5, 5)) # y座標軸的刻度範圍
24     # plt.xlabel("預測偏差數: %d" % (n_errors))
25     legend = plt.legend(loc='upper left')
26     legend.legendHandles[0]._sizes = [10]
27     plt.show()

　　create_train()建立了一個數據集，這個數據集中包含了兩簇正常的樣本和20個可能的異常樣本，之因此是「可能的異常樣本」，是因爲X_outliers中的數據是隨機採樣的，可能靠近正常的樣本簇：

　　接下來按照前面的介紹依次實現數據集中每個樣本的k距離、k距離鄰域、可達距離、局部可達密度和局部異常因子，並修改display()，將每一個樣本點的異常因子標記在圖像上。完整代碼：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 
 4 def create_train():
 5     np.random.seed(42) # 設置seed使每次生成的隨機數都相等
 6     # 生成100個2維數據，它們是以0爲均值、以1爲標準差的正態分佈
 7     X_inliers = 0.3 * np.random.randn(100, 2)
 8     # 構造兩組間隔必定距離的樣本點做爲訓練數據
 9     X_inliers = np.r_[X_inliers + 2, X_inliers - 2]
10     # 構造20個可能的異常數據，從一個均勻分佈[low,high)中隨機採樣
11     X_outliers = np.random.uniform(low=-4, high=4, size=(20, 2))
12     # 將X_inliers和X_outliers鏈接起來做爲訓練集
13     return np.r_[X_inliers, X_outliers]
14 
15 def display(X):
16     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
17     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
18     plt.title("局部異常因子(LOF)")
19     plt.scatter(X[:, 0], X[:, 1], color='k', s=3., label='樣本點')
20 
21     max_lof, min_lof = max(LOF), min(LOF)
22     # 每一個樣本點的半徑，樣本點越異常，半徑越大(歸一化處理)
23     radius = [1000 * ((x - min_lof) / (max_lof - min_lof)) for x in LOF]
24     plt.scatter(X[:, 0], X[:, 1], s=radius, edgecolors='r', facecolors='none', label='異常度')
25 
26     plt.axis('tight') # 讓座標軸的比例尺適應數據量
27     plt.xlim((-5, 5)) # x座標軸的刻度範圍
28     plt.ylim((-5, 5)) # y座標軸的刻度範圍
29     # plt.xlabel("預測偏差數: %d" % (n_errors))
30     legend = plt.legend(loc='upper left')
31     legend.legendHandles[0]._sizes = [10]
32     legend.legendHandles[1]._sizes = [20]
33     plt.show()
34 
35 D_K = [] # X中每一個樣本點的k距離
36 NEIGHBORS = [] # X中每一個樣本點的k距離鄰域
37 LRD = [] # X中每一個樣本點的局部可達密度
38 
39 def euclid(x_i, x_j):
40     # 兩點間的歐幾里得距離
41     return np.linalg.norm(x_i - x_j)
42 
43 def all_D_k(X, k):
44     for x_i in X:
45         # x_i距全部樣本點的歐幾里得距離
46         x_and_d = [(x_j, euclid(x_i, x_j)) for x_j in X if (x_i != x_j).any()]
47         x_and_d = sorted(x_and_d, key=lambda x: x[1])  # 按照距離排序
48         D_K.append(x_and_d[k - 1][1]) # x_i的k距離，此處假設x_and_d中沒有重複的距離
49         NEIGHBORS.append([x for x, d in x_and_d[:k]])  # x_i的k距離鄰域
50     for x_i in X: # 計算每一個樣本點的局部可達密度
51         neighbors = neighbors_k(x_i, X, k)  # x_i的k鄰域
52         LRD.append(1 / (sum([rd(x_t, x_i, X, k) for x_t in neighbors]) / len(neighbors)))
53 
54 def d_k(x, X, k):
55     # x的k距離
56     i = np.argwhere(X == x)[0, 0]  # x在訓練樣本中的索引
57     return D_K[i]
58 
59 def neighbors_k(x, X, k):
60     # x的k距離鄰域
61     i = np.argwhere(X == x)[0,0] # x在訓練樣本中的索引
62     return NEIGHBORS[i]
63 
64 def rd(x_i, x_j, X, k):
65     # x_i到x_j的可達距離
66     return max(d_k(x_i, X, k), euclid(x_i, x_j))
67 
68 def lrd(x, X, k):
69     # x的局部可達密度
70     i = np.argwhere(X == x)[0, 0]  # x在訓練樣本中的索引
71     return LRD[i]
72 
73 def lof(X, k):
74     # 數據集中每一個數據樣本的局部異常因子
75     LOF = []
76     for i, x_i in enumerate(X):
77         neighbors = neighbors_k(x_i, X, k)  # x_i的k鄰域
78         LOF.append(sum([lrd(x_t, X, k) for x_t in neighbors]) / len(neighbors) / lrd(x_i, X, k))
79     return LOF
80 
81 if __name__ == '__main__':
82     X = create_train()
83     k = 20
84     all_D_k(X, k) # 計算X中每一個樣本點的第k近的距離
85     LOF = lof(X,k)
86     display(X)
87 
88     for i, lof in enumerate(LOF):
89         print(lof, '\t', end='')
90         if (i + 1) % 7 == 0:
91             print()