異常檢測(2)——基於機率統計的異常檢測（1）

　　某個工廠生產了一批手機屏幕，爲了評判手機屏幕的質量是否達到標準，質檢員須要收集每一個樣本的若干項指標，好比大小、質量、光澤度等，根據這些指標進行打分，最後判斷是否合格。如今爲了提升效率，工廠決定使用智能檢測進行第一步篩選，質檢員只須要重點檢測被系統斷定爲「不合格」的樣本。

　　智能檢測程序須要根據大量樣本訓練一個函數模型，也許咱們的第一個想法是像監督學習那樣，爲樣本打上「正常」和「異常」的標籤，而後經過分類算法訓練模型。假設x_test是數據樣本，predict(x_test)來判斷x_test是不是合格樣本。某個偷懶的傢伙寫下了這樣的代碼：

 def predict(xtest):
        return 1

　　因爲工廠的質量管理過硬，僅有極少數不合格樣本，所以這段荒唐的預測竟然展示出極高的準確率！這是因爲嚴重的數據偏斜致使的，或許咱們能夠經過查準率（Precision）和召回率（Recall）兩個指標識別出這段不負責任的代碼，可是當你再次試圖使用某個監督學習算法時，仍然會面對一樣的問題——僅有極少數不合格樣本，以致於監督學習沒法學到足夠的知識。可否從極度偏斜的數據中學習出一個有效的檢測模型呢？固然能，這就是基於統計的異常檢測。這類方法一般會假設給定的數據集服從一個隨機分佈模型，將與模型不一致的樣本視爲異常樣本。其中最經常使用的兩種分佈模型是一元正態分佈模型和多元正態分佈模型。

算法模型

　　在正態分佈的假設下，若是有一個新樣本x，當x的正態分佈值小於某個閾值時，就能夠認爲這個樣本是異常的。

　　在正態分佈中，μ-3σ<=x<=μ+3σ的區域包含了絕大部分數據，能夠以此爲參考，調整ε的值：

　　如今有一個包含了m個一維數據的訓練集：

　　能夠經過下面的函數判斷一個樣本是不是異常的：

　　這裏x⁽ⁱ⁾是已知的，μ和σ纔是未知的，咱們的目的是設法根據訓練集求得μ和σ的值，以獲得一個肯定的函數模型。具體來講，經過最大似然估計量能夠得出下面的結果：

　　具體推導過程參考 機率筆記11——一維正態分佈的最大似然估計

 

算法實現

　　咱們經過一些模擬數據來一窺異常檢測算法的究竟。

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 
  4 def create_data():
  5     '''
  6     建立訓練數據和測試數據
  7     :return: X_train:訓練集， X_test:測試集
  8     '''
  9     np.random.seed(42)  # 設置seed使每次生成的隨機數都相等
 10     m, s = 3, 0.1 # 設置均值和方差
 11     X_train = np.random.normal(m, s, 100) # 100個一元正態分佈數據
 12     # 構造10測試數據，從一個均勻分佈[low,high)中隨機採樣
 13     X_test = np.random.uniform(low=m - 1, high=m + 1, size=10)
 14     return X_train, X_test
 15 
 16 def plot_data(X_train, X_test):
 17     '''
 18     數據可視化
 19     :param X_train: 訓練集
 20     :param X_test: 測試集
 21     :return:
 22     '''
 23     fig = plt.figure(figsize=(10, 4))
 24     plt.subplots_adjust(wspace=0.5)  # 調整子圖之間的左右邊距
 25     fig.add_subplot(1, 2, 1)  # 繪製訓練數據的分佈
 26     plt.scatter(X_train, [0] * len(X_train), color='blue', marker='x', label='訓練數據')
 27     plt.title('訓練數據的分佈狀況')
 28     plt.xlabel('x')
 29     plt.ylabel('y')
 30     plt.legend(loc='upper left')
 31 
 32     fig.add_subplot(1, 2, 2)  # 繪製總體數據的分佈
 33     plt.scatter(X_train, [0] * len(X_train), color='blue', marker='x', label='訓練數據')
 34     plt.scatter(X_test, [0] * len(X_test), color='red', marker='^',label='測試數據')
 35     plt.title('總體數據的分佈狀況')
 36     plt.xlabel('x')
 37     plt.ylabel('y')
 38     plt.legend(loc='upper left')
 39 
 40     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
 41     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
 42     plt.show()
 43 
 44 def fit(X_train):
 45     '''
 46     擬合數據，訓練模型
 47     :param X_train: 訓練集
 48     :return:  mu:均值, sigma:方差
 49     '''
 50     global mu, sigma
 51     mu = np.mean(X_train)  # 計算均值μ
 52     sigma = np.var(X_train) # 計算方差 σ^2
 53 
 54 def gaussian(X):
 55     '''
 56     計算正態分佈
 57     :param X: 數據集
 58     :return: 數據集的密度值
 59     '''
 60     return np.exp(-((X - mu) ** 2) / (2 * sigma)) / (np.sqrt(2 * np.pi) * np.sqrt(sigma))
 61 
 62 def get_epsilon(n=3):
 63     ''' 調整ε的值，默認ε=3σ '''
 64     return np.sqrt(sigma) * n
 65 
 66 def predict(X):
 67     '''
 68     檢測訓練集中的數據是不是正常數據
 69     :param X: 待預測的數據
 70     :return: P1:數據的密度值, P2:數據的異常檢測結果，True正常，False異常
 71     '''
 72     P1 = gaussian(X) # 數據的密度值
 73     epsilon = get_epsilon()
 74     P2 = [p > epsilon for p in P1] # 數據的異常檢測結果，True正常，False異常
 75     return P1, P2
 76 
 77 def plot_predict(X):
 78     '''可視化異常檢測結果 '''
 79     epsilon = get_epsilon()
 80     xs = np.linspace(mu - epsilon, mu + epsilon, 50)
 81     ys = gaussian(xs)
 82     plt.plot(xs, ys, c='g', label='擬合曲線')  # 繪製正態分佈曲線
 83 
 84     P1, P2 = predict(X)
 85     normals_idx = [i for i, t in enumerate(P2) if t == True] # 正常數據的索引
 86     plt.scatter([X[i] for i in normals_idx], [P1[i] for i in normals_idx],
 87                 color='blue', marker='x', label='正常數據')
 88     outliers_idx = [i for i, t in enumerate(P2) if t == False] # 異常數據的索引
 89     plt.scatter([X[i] for i in outliers_idx], [P1[i] for i in outliers_idx],
 90                 color='red', marker='^', label='異常數據')
 91     plt.title('檢測結果，共有{}個異常數據'.format(len(outliers_idx)))
 92     plt.xlabel('x')
 93     plt.ylabel('y')
 94     plt.legend(loc='upper left')
 95     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
 96     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
 97     plt.show()
 98 
 99 if __name__ == '__main__':
100     mu, sigma = 0, 0 # 模型的均值μ和方差σ^2
101     X_train, X_test = create_data()
102     plot_data(X_train, X_test)
103     fit(X_train)
104     print('μ = {}, σ^2 = {}'.format(mu, sigma))
105     plot_predict(np.r_[X_train, X_test])

　　create_data()建立了100個知足X~N(μ，σ²)的數據做爲訓練集；10個在 之間均勻分佈的數據做爲測試數據，它們極有多是異常數據。plot_data()可視化了訓練集，再把測試數據加進去一併展現：

　　能夠看出，大部分訓練數據都集中在正態分佈的均值區域，而異常數據偏向於「倒鍾」的兩端。

　　接下來使用fit()方法對異常檢測模型進行訓練，獲得的結果是μ = 2.98961534826059, σ^2 = 0.008165221946938589。

　　獲得了模型參數後就可使用目標函數對數據進行預測。gussian(X)實現了正態分佈的密度函數；predict(X)將對X中的全部樣本進行檢測，並返回X對應的檢測結果列表。其可視化結果是：

一元模型的問題

　　在面對多維數據時，基於一元正態分佈的異常檢測能夠單獨抽取某一維度進行檢測，一般也能工做的很好，但這裏有一個假設—全部維度都符合正態分佈，而且各維度都是獨立的，若是兩個維度之間存在相關性，那麼基於一元正態分佈的異常檢測就可能會出現很大程度的誤判。

　　人的身高和鞋碼存在着關聯關係，通常來講，身高是腳長的7倍左右。假設某地區成年男子的身高符合μ=1.7，σ²=0.036的正態分佈，咱們用下面的代碼模擬身高和腳長的數據。

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 
 4 def create_train():
 5     '''
 6     構造2維訓練集
 7     :return: X1:第1緯度的列, X2:第2維度的列表
 8     '''
 9     np.random.seed(21)  # 設置seed使每次生成的隨機數都相等
10     mu, sigma = 1.7, 0.036  # 設置均值和方差
11     X1 = np.random.normal(mu, sigma, 200)  # 生成200個符合正態分佈的身高數據
12     # 設置身高對應的鞋碼，鞋碼=身高/7±0.02
13     X2 = (X1 / 7) + np.random.uniform(low=-0.01, high=0.01, size=len(X1))
14     return X1, X2
15 
16 def plot_train(X1, X2):
17     '''
18     可視化訓練集
19     :param X1: 訓練集的第1維度
20     :param X2: 訓練集的第2維度
21     '''
22     fig = plt.figure(figsize=(10, 4))
23     plt.subplots_adjust(hspace=0.5)  # 調整子圖之間的上下邊距
24     # 數據的散點圖
25     fig.add_subplot(2, 1, 1)
26     plt.scatter(X1, X2, color='k', s=3., label='訓練數據')
27     plt.legend(loc='upper left')
28     plt.xlabel('身高')
29     plt.ylabel('腳長')
30     plt.title('數據分佈')
31     # x1維度的直方圖
32     fig.add_subplot(2, 2, 3)
33     plt.hist(X1, bins=40)
34     plt.xlabel('身高')
35     plt.ylabel('頻度')
36     plt.title('身高直方圖')
37     # x2維度的直方圖
38     fig.add_subplot(2, 2, 4)
39     plt.hist(X2, bins=40)
40     plt.xlabel('腳長')
41     plt.ylabel('頻度')
42     plt.title('腳長直方圖')
43 
44     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
45     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
46     plt.show()
47 
48 
49 if __name__ == '__main__':
50     X1_train, X2_train = create_train()
51     plot_train(X1_train, X2_train)

　　create_train()方法先是構造了符合正態分佈的身高數據，以後根據身高構造腳長，併爲腳長加入了在±0.01之間的噪聲，以便使腳長產生一些震盪。plot_train(X1, X2)用散點圖和直方圖的形式將數據可視化：

　　從上面的直方圖能夠看出，身高和與之存在線性相關性的腳長都符合正態分佈。咱們能夠用前面的代碼的計算兩個維度的均值和方差，在根據這些參數畫出每一個維度的擬合曲線。

　　就訓練集而言，低於1.6米或高於1.8米的身高被認爲偏向於異常值，與這些身高有關聯關係的腳長的數值也被認爲偏向於異常值。如今有一個1.73米、腳長是0.23米的人，這個小腳男士遠離了大部分樣本點，應當被視爲異常數據：

　　如今的問題是，對於身高的分佈來講，1.73米是一個大衆化的身高，對於腳長來講，0.23米也不算太差：

　　不管從身高維度仍是腳長維度，這個小腳男士都被認爲是正常數據，這與咱們指望的結果不符。其緣由是一元正態分佈分別從兩個獨立的維度看待問題，而忽略了這兩個維度間的相關性。爲了識別相關性，咱們須要使用更高維度的正態分佈模型。（下篇繼續）

---

　　做者：我是8位的

　　出處：http://www.cnblogs.com/bigmonkey

　　本文以學習、研究和分享爲主，如需轉載，請聯繫本人，標明做者和出處，非商業用途！ 

　　掃描二維碼關注公做者衆號「我是8位的」