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

　　書接上文，繼續討論基於多元正態分佈的異常檢測算法。

　　

　　如今有一個包含了m個數據的訓練集，其中的每一個樣本都是一個n維數據：

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

　　咱們的目的是設法根據訓練集求得μ和σ，以獲得一個肯定的多元分正態布模型。具體來講，經過最大似然估計量能夠得出下面的結論：

　　其中Σ是協方差對角矩陣，最終求得的多元正態分佈模型能夠寫成：

　　關於最大似然估計量、協方差矩陣和多元正態分佈最大似然估計的具體推導過程，可參考：

　　機率筆記12——多維正態分佈的最大似然估計

　　機率筆記10——矩估計和最大似然

 

　

　　代碼：

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 from mpl_toolkits.mplot3d import Axes3D
  4
  5 def create_data(model='train', count=200):
  6     '''
  7     構造2維訓練集
  8     :param model: train:訓練集,  test:測試集
  9     :param count: 樣本數量
 10     :return: X1:第1緯度的列, X2:第2維度的列表
 11     '''
 12     np.random.seed(21)  # 設置seed使每次生成的隨機數都相等
 13     X1 = np.random.normal(1.7, 0.036, count)  # 生成200個符合正態分佈的身高數據
 14     low, high = -0.01, 0.01
 15     if model == 'test':
 16         low, high = -0.05, 0.05
 17     # 設置身高對應的鞋碼，正常鞋碼=身高/7±0.01
 18     X2 = X1 / 7 + np.random.uniform(low=low, high=high, size=len(X1))
 19     return X1, X2
 20
 21 def plot_train(X1, X2):
 22     '''
 23     可視化訓練集
 24     :param X1: 訓練集的第1維度
 25     :param X2: 訓練集的第2維度
 26     '''
 27     fig = plt.figure(figsize=(10, 4))
 28     plt.subplots_adjust(hspace=0.5)  # 調整子圖之間的上下邊距
 29     # 數據的散點圖
 30     fig.add_subplot(2, 1, 1)
 31     plt.scatter(X1, X2, color='k', s=3., label='訓練數據')
 32     plt.legend(loc='upper left')
 33     plt.xlabel('身高')
 34     plt.ylabel('腳長')
 35     plt.title('數據分佈')
 36     # 身高維度的直方圖
 37     fig.add_subplot(2, 2, 3)
 38     plt.hist(X1, bins=40)
 39     plt.xlabel('身高')
 40     plt.ylabel('頻度')
 41     plt.title('身高直方圖')
 42     # 腳長維度的直方圖
 43     fig.add_subplot(2, 2, 4)
 44     plt.hist(X2, bins=40)
 45     plt.xlabel('腳長')
 46     plt.ylabel('頻度')
 47     plt.title('腳長直方圖')
 48
 49     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
 50     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
 51     plt.show()
 52
 53 def fit(X_train):
 54     '''
 55     擬合數據，訓練模型
 56     :param X_train: 訓練集
 57     :return:  mu:均值, sigma:方差
 58     '''
 59     global mu, sigma
 60     X = np.mat(X_train.T)
 61     m, n = X.shape
 62     mu = np.mean(X, axis=1)  # 計算均值μ,axis=1表示對每個子數組計算均值
 63     sigma = np.mat(np.cov(X)) # 計算Σ,等同於(X - mu) * (X - mu).T / len(X_train)
 64
 65 def gaussian(X_test):
 66     '''
 67     計算正態分佈
 68     :param X_test: 測試集
 69     :return: 數據集的密度值
 70     '''
 71     global mu, sigma
 72     m, n = np.shape(X_test)
 73     sig_det = np.linalg.det(sigma)  # 計算det(Σ)
 74     sig_inv = np.linalg.inv(sigma)  # Σ的逆矩陣
 75     r = []
 76     for x in X_test:
 77         x = np.mat(x).T - mu
 78         g = np.exp(-x.T * sig_inv * x / 2) *  ((2 * np.pi) ** (-n / 2) * (sig_det ** (-0.5)))
 79         r.append(g[0, 0])
 80     return r
 81
 82 def sel_epsilon(X_train):
 83     '''
 84     選擇合適的ε
 85     :param X_train:
 86     :return: ε
 87     '''
 88     g_val = gaussian(X_train)
 89     return np.min(g_val) + 0.0001
 90
 91 def predict(X):
 92     '''
 93     檢測訓練集中的數據是不是正常數據
 94     :param X: 待預測的數據
 95     :return: P1:數據的密度值, P2:數據的異常檢測結果，True正常，False異常
 96     '''
 97     P1 = gaussian(X) # 數據的密度值
 98     P2 = [p > epsilon for p in P1] # 數據的異常檢測結果，True正常，False異常
 99     return P1, P2
100
101 def plot_predict(X):
102     '''可視化異常檢測結果 '''
103     P1, P2 = predict(X)
104     normals_idx = [i for i, t in enumerate(P2) if t == True]  # 正常數據的索引
105     outliers_idx = [i for i, t in enumerate(P2) if t == False]  # 異常數據的索引
106     normals_x = np.array([X[i] for i in normals_idx])  # 正常數據
107     outliers_x =  np.array([X[i] for i in outliers_idx])  # 異常數據
108
109     fig1 = plt.figure(num='fig1') # 散點圖
110     ax = Axes3D(fig1)
111     ax.scatter(normals_x[:,0], normals_x[:,1],
112                [P1[i] for i in normals_idx], label='正常數據')
113     ax.scatter(outliers_x[:,0], outliers_x[:,1],
114                [P1[i] for i in outliers_idx], c='r', marker='^', label='異常數據')
115     ax.set_title('共有{}個異常數據'.format(len(outliers_idx)))
116     ax.axis('tight')  # 讓座標軸的比例尺適應數據量
117     ax.set_xlabel('身高')
118     ax.set_ylabel('腳長')
119     ax.set_zlabel('p(x)')
120     ax.legend(loc='upper left')
121
122     n = 100
123     xs, ys = np.meshgrid(np.linspace(min(X1_train), max(X1_train), n),
124                          np.linspace(min(X2_train), max(X2_train), n))
125     zs = [gaussian(np.c_[xs[i], ys[i]]) for i in range(n)]
126
127     fig2 = plt.figure(num='fig2')
128     ax = Axes3D(fig2)
129     # 3維空間的擬合曲面
130     ax.plot_surface(xs, ys, zs, alpha=0.5, cmap=plt.get_cmap('rainbow'))
131     ax.scatter(normals_x[:, 0], normals_x[:, 1],
132                [P1[i] for i in normals_idx], label='正常數據')
133     ax.scatter(outliers_x[:, 0], outliers_x[:, 1],
134                [P1[i] for i in outliers_idx], c='r', marker='^', label='異常數據')
135     ax.axis('tight')  # 讓座標軸的比例尺適應數據量
136     ax.set_xlabel('身高')
137     ax.set_ylabel('腳長')
138     ax.set_zlabel('p(x)')
139     ax.legend(loc='upper left')
140
141     fig3 = plt.figure(num='fig3')
142     plt.scatter(normals_x[:, 0], normals_x[:, 1], s=30., c='k', label='正常數據')
143     plt.scatter(outliers_x[:, 0], outliers_x[:, 1], c='r', marker='^', s=30., label='異常數據')
144     plt.contour(xs, ys, zs, 80, alpha=0.5) # 等高線
145     plt.axis('tight')  # 讓座標軸的比例尺適應數據量
146     plt.xlabel('身高')
147     plt.ylabel('腳長')
148     plt.legend(loc='upper left')
149
150     plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
151     plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
152     plt.show()
153
154 if __name__ == '__main__':
155     X1_train, X2_train = create_data()
156     plot_train(X1_train, X2_train)
157
158     X_train = np.c_[X1_train, X2_train]
159     mu, sigma = np.mat([]), np.mat([])
160     fit(X_train)
161     epsilon = sel_epsilon(X_train)
162
163     X_test = np.c_[create_data(model='test', count=20)]
164     X = np.r_[X_train, X_test]
165     plot_predict(X)

　　爲了簡單起見，咱們認爲X_train 中的數據所有是正常數據。fit(X_train)計算多元正態分佈的模型參數，gaussian(X_test)根據目標函數計算樣本的多元正態分佈密度值。在知道了算法原理的請看下，fit(X_train)和gaussian(X_test)都毫無神祕性可言。

　　predict(X)將對X中的全部樣本進行檢測，並返回X對應的檢測結果列表，列表中的元素是一個二元組，第一個元素記錄x⁽ⁱ⁾是不是正常數據，第二個元素記錄p(x⁽ⁱ⁾;μ,Σ)。因爲已經假設了X_train中所有是正常數據，所以這裏選擇X_train中最小的密度值做爲ε。

　　X_test中的20個測試數據是可能的異常樣本。plot_predict(X)展現了空間樣本點、空間擬合曲面和等高線：

　　

　　

---

　　做者：我是8位的

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

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

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