異常檢測(Anomaly detection)是機器學習中一種常見應用，其目標是識別資料集中的異常或不尋常模式。儘管通常被歸類為非監督學習問題，異常檢測卻具有與監督學習相似的特徵。在異常檢測中，我們通常處理的是未標記的資料，即沒有明確的標籤指示哪些樣本是異常的。相反，演算法需要根據資料本身的特徵來確定異常。這使得異常檢測成為一項挑戰，因為異常通常是稀有事件，不易獲取大量標記的異常資料以進行訓練。儘管異常檢測被歸類為非監督學習，但我們可以將其視為一種“半監督”學習，因為我們通常有一些正常樣本，但沒有足夠的異常樣本。這種情況下，我們可以利用正常樣本來構建模型，然後將其應用於整個資料集以檢測異常。
一、問題來源
假想你是一個飛機引擎製造商，當你生產的飛機引擎從生產線上流出時，你需要進行QA(質量控制測試)，而作為這個測試的一部分，你測量了飛機引擎的一些特徵變數，比如引擎運轉時產生的熱量，或者引擎的振動等。如果你生產了\(m\)個引擎的話，那麼你就有了一個資料集，從\(X^{(1)}\)到\(X^{(m)}\)，將這些資料繪製成圖表，如下圖1中紅色X所示。這裡每個叉，都是無標籤的資料。那麼，異常檢測問題可以定義為：我們假設後來有一天，你有一個新的飛機引擎從生產線上流出，而你的新飛機引擎有特徵變數\(x_{test}\)，我們希望知道這個新的飛機引擎是否有某種異常，或者說，我們希望判斷這個引擎是否需要進一步測試。因為，如果它看起來像一個正常的引擎，那麼我們可以直接將它運送到客戶那裡，而不需要進一步的測試。從直覺上先來理解下異常檢測問題，現在假設，這個新的飛機引擎，如下圖2，它出現在比較靠近中心的位置（綠X），那麼它看起來是 OK 的，但如果它出現的位置距離資料集的中心比較偏遠，那麼它有可能是 anomaly，我們可能需要對它再進行測試一下。
| 圖1 | 圖2 |
|---|---|
 |
 |
讓我們定義這個問題：
假設我們有資料集 \(x^{(1)},x^{(2)},..,x^{(m)}\)，通常我們假定 \(m\) 個樣本都是正常的，然後我們希望有一個演算法告訴我們，一個新的樣本資料 \(x_{test}\) 是否異常，我們要採取的做法是，給定無標籤的訓練集，我們將對資料建模，即 \(p(x)\)，也就是說，我們將對 \(x\)的機率分佈建模，其中 這些\(x\) 是特徵變數，如引擎運轉時產生的熱量，或者引擎的振動等。因此，當我們建立了 \(x\) 的機率模型之後，我們就有對於一個新的飛機引擎 \(x_{test}\)，它的機率\(p(x_{test})\)如果低於一個閾值 \(\varepsilon\)，那麼就將其標記為異常（anomaly）。這意味著，對於整個資料集，這個點出現的機率非常的低。反之，如果機率大於等於給定的閾值 \(\varepsilon\)，我們就認為它是正常的（normal）。這種方法稱為密度估計，可表示為：$$\quad p(x_{test}) \begin{cases} < \varepsilon & Anomaly \\ \geq \varepsilon & Normal \end{cases}$$
如下圖3，給定圖中這樣的訓練集，如果你建立了一個模型，你將很可能發現飛機引擎，即模型\(p(x)\) 會認為在中心區域的這些點機率相當大，而稍微遠離中心區域的點機率會小點，更遠地方的點，機率會更小，而在藍色圈外面的點，將成為異常點。
| 圖3 | 圖4 |
|---|---|
 |
 |
異常檢測演算法的實際應用有很多，最常見的應用是欺詐檢測。假設有很多使用者，每個使用者都在從事不同的活動，也許是在個人網站上，也許是在一個實體工廠之類的地方，可以透過對不同的使用者活動計算特徵變數，如\(x_1\)是使用者登陸的頻率，\(x_2\)也許是使用者訪問某個頁面的次數或者交易次數，\(x_3\)是使用者在論壇上發貼的次數，\(x_4\)是使用者的打字速度等等，然後建立一個模型\(p(x)\)，這時可以用它來發現網站上的行為異常的使用者。只需要看哪些使用者的\(p(x)\)機率小於\(\varepsilon\)，接下來拿來這些使用者的檔案做進一步篩選，或者要求這些使用者驗證他們的身份等等，從而讓你的網站防禦異常行為或者欺詐盜號等行為。
異常檢測的另一個例子是在工業生產領域，事實上我們之前所說的飛機引擎的問題，透過模型\(p(x)\)找到異常的飛機引擎，然後要求進一步細查這些引擎的質量。 此外對於資料中心的計算機監控也是異常檢測的應之一，如果你管理一個計算機叢集或者一個資料中心，其中有許多計算機，我們可以為每臺計算機計算特徵變數如：計算機的記憶體消耗、硬碟訪問量、CPU負載或者一些更加複雜的特徵，例如一臺計算機的CPU負載與網路流量的比值，那麼給定正常情況下資料中心中計算機的特徵變數，可以建立\(p(x)\)模型，如果有一個計算機它的機率\(p(x)\)非常小，這時可以認為這個計算機執行不正常，從而進一步要求系統管理員檢視其工作狀況，目前這種技術實際正在被各大資料中心使用用來監測大量計算機可能發生的異常。
二、基於統計學的異常檢測
異常點（outlier）是一個資料物件，它明顯不同於其他的資料物件，就好像它是被不同的機制產生的一樣。
2.1 正態分佈一元離群點的檢測方法
一元正態分佈函式：$$p(x) = {1 \over {\sqrt {2\pi } \sigma }}{e^{ - {{{{(x - u)}^2}} \over {2{\sigma ^2}}}}}$$
其中，\(\mu\)為資料的均值，\(\sigma\)為資料的標準差，\(\sigma\)越小，對應的影像越尖。
假設有\(n\) 個點 \((x_{1},...,x_{n})\)，那麼可以估計出這\(n\)個點的均值\(\mu\) 和方差\(\sigma\)。均值和方差分別被定義為：
在正態分佈的假設下，區域 \(\mu\pm 3\sigma\)包含了99.7% 的資料，如果某個值距離分佈的均值\(\mu\)超過了\(3\sigma\)，那麼這個值就可以被簡單的標記為一個異常點（outlier）。
2.2 多元正態分佈離群點的檢測方法
涉及兩個或者兩個以上變數的資料稱為多後設資料，很多一元離群點的檢測方法都可以擴充套件到高維空間中，從而處理多後設資料。
- 基於一元正態分佈的離群點檢測方法
假設\(n\)維的資料集合形如 \(x_{i}=(x_{i,1},...,x_{i,n}), i\in \{1,...,m\}\)，那麼可以計算每個維度的均值和方差\(\mu_{j},\sigma_{j}, j\in\{1,...,n\}\)。 具體來說，對於j\in {1,...,n}，可以計算
在正態分佈的假設下，如果有一個新的資料 \(x\)，可以計算機率\(p(x)\) 如下：
根據機率值的大小就可以判斷 \(x\) 是否屬於異常值。
- 多維正態分佈——高斯分佈的離群點檢測方法
假設 \(n\) 維的資料集合 \(x=(x_{1},...,x_{n})\)，可以計算\(n\)維的均值向量
和 \(n\times n\) 的協方差矩陣：
如果有一個新的資料 \(x_{n+1}=(x_{n+1,1},...,x_{n+1,n})\)，可以計算我們首先計算所有特徵的平均值，然後再估計協方差矩陣：
注意\(\mu\)是一個向量，其每一個單元都是原特徵矩陣中一行資料的均值。最後我們計算多元高斯分佈的$ p(x)$：
其中：\(|\Sigma|\) 是 \(\Sigma\) 的行列式；$ \Sigma^{-1}$ 是 $\Sigma $的逆矩陣。
異常判定：當 $p(x) < \epsilon $ 時，判定為異常，其中$ \epsilon $ 是預先設定的閾值。如果這個分數低於某個閾值 $ \epsilon $，那麼該樣本就被認為是異常的。
import numpy as np
import matplotlib.pyplot as plt
# 生成虛構的二維資料
np.random.seed(0)
mean = [0, 0]
cov = [[1, 0.5], [0.5, 1]]
x = np.random.multivariate_normal(mean, cov, 300)
# 新資料點
new_point = np.array([7, 6])
# 視覺化生成的資料及新資料點
plt.figure(figsize=(8, 5))
plt.scatter(x[:, 0], x[:, 1], edgecolors='b', label='Data')
plt.scatter(new_point[0], new_point[1], color='r', marker='x', label='New Point (7, 6)')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('Generated Data with New Point')
plt.legend()
plt.show()
def estimate_gaussian(x):
"""估計高斯分佈的均值和方差"""
mu = x.mean(axis=0) # 求每列的均值
sigma2 = x.var(axis=0) # 求每列的方差，這裡自由度為m
return mu, sigma2
def multivariate_gaussian(x, mu, sigma2):
"""多元高斯分佈密度函式"""
p = np.zeros((x.shape[0], 1))
n = len(mu)
if np.ndim(sigma2) == 1:
sigma2 = np.diag(sigma2) # 對角陣
for i in range(x.shape[0]):
p[i] = (2 * np.pi) ** (-n / 2) * np.linalg.det(sigma2) ** (-1 / 2) * np.exp(
-0.5 * (x[i, :] - mu).T @ np.linalg.inv(sigma2) @ (x[i, :] - mu))
return p
# 估計高斯分佈的引數
mu, sigma2 = estimate_gaussian(x)
# 計算機率密度
p = multivariate_gaussian(x, mu, sigma2)
# 計算新資料點的機率密度
p_new_point = multivariate_gaussian(new_point.reshape(1, -1), mu, sigma2)
# 選擇最優閾值
def select_threshold(p):
"""選擇最優閾值"""
bestF1 = 0
bestEpsilon = 0
for epsilon in np.linspace(min(p), max(p), 1001):
y_predict = np.zeros(p.shape)
y_predict[p < epsilon] = 1 # 把小於閾值的設為1(異常)
tp = np.sum(y_predict) # 真正例
precision = tp / np.sum(y_predict == 1)
recall = tp / len(p)
F1 = (2 * precision * recall) / (precision + recall)
if F1 > bestF1:
bestF1 = F1
bestEpsilon = epsilon
return bestEpsilon, bestF1
epsilon, F1 = select_threshold(p)
print("最優閾值:", epsilon)
print("最優F1分數:", F1)
# 解釋最優F1分數
print("\nF1 分數是精確率（Precision）和召回率（Recall）的調和平均。")
print("精確率衡量了模型在預測為異常時的準確性，即異常點中真正的異常點的比例。")
print("召回率衡量了模型能夠發現真正的異常點的能力，即所有真正異常點中被發現的比例。")
print("F1 分數越高，說明模型在精確率和召回率之間取得了更好的平衡。")
print("在這個例子中，最優 F1 分數為", F1, "，表示模型在判定異常點時達到了較好的平衡。")
# 判斷新資料點是否為異常點
if p_new_point < epsilon:
print("新資料點 (7, 6) 是異常點。")
else:
print("新資料點 (7, 6) 不是異常點。")
2.3 \(\chi^{2}\) 統計量檢測多元離群點
在正態分佈的假設下，\(\chi^{2}\) 統計量可以用來檢測多元離群點。對於某個物件\(\bold{a}\)，\(\chi^{2}\) 統計量是
其中，\(a_{i}\)是\(\bold{a}\)在第 \(i\) 維上的取值，\(E_{i}\)是所有物件在第 \(i\) 維的均值，\(n\)是維度。如果物件\(\bold{a}\)的\(\chi^{2}\) 統計量很大，那麼該物件就可以認為是離群點。
三、其他檢查方法
3.1 基於聚類的方法
基於聚類的異常檢測方法通常依賴下列假設，1）正常資料例項屬於資料中的一個簇，而異常資料例項不屬於任何簇；2）正常資料例項靠近它們最近的簇質心，而異常資料離它們最近的簇質心很遠；3）正常資料例項屬於大而密集的簇，而異常資料例項要麼屬於小簇，要麼屬於稀疏簇；透過將資料歸分到不同的簇中，異常資料則是那些屬於小簇或者不屬於任何一簇或者遠離簇中心的資料。
將距離簇中心較遠的資料作為異常點：這類方法有 SOM、K-means、最大期望( expectation maximization，EM)及基於語義異常因子( semantic anomaly factor)演算法等；將聚類所得小簇資料作為異常點：代表方法有K-means聚類；將不屬於任何一簇作為異常點：代表方法有 DBSCAN、ROCK、SNN 聚類。
3.2 基於深度的方法
該方法將資料對映到 \(k\) 維空間的分層結構中，並假設異常值分佈在外圍，而正常資料點靠近分層結構的中心（深度越高）。半空間深度法( ISODEPTH 法) ，透過計算每個點的深度，並根據深度值判斷異常資料點。最小橢球估計 ( minimum volume ellipsoid estimator，MVE)法。根據大多數資料點( 通常為 ＞ 50% ) 的機率分佈模型擬合出一個實線橢圓形所示的最小橢圓形球體的邊界，不在此邊界範圍內的資料點將被判斷為異常點。
| 圖5 | 圖6 |
|---|---|
 |
 |
3.3 基於神經網路的方法
代表方法有自動編碼器( autoencoder，AE) ，長短期記憶神經網路（LSTM）等。LSTM可用於時間序列資料的異常檢測：利用歷史序列資料訓練模型，檢測與預測值差異較大的異常點。Autoencoder異常檢測 Autoencoder本質上使用了一個神經網路來產生一個高維輸入的低維表示。Autoencoder與主成分分析PCA類似，但是Autoencoder在使用非線性啟用函式時克服了PCA線性的限制。演算法的基本上假設是異常點服從不同的分佈。根據正常資料訓練出來的Autoencoder，能夠將正常樣本重建還原，但是卻無法將異於正常分佈的資料點較好地還原，導致其基於重構誤差較大。當重構誤差大於某個閾值時，將其標記為異常值。
四、信用卡反欺詐
專案為kaggle上經典的信用卡欺詐檢測，該資料集質量高，正負樣本比例非常懸殊。我們在此專案主要用了無監督的Autoencoder新穎點檢測，根據重構誤差識別異常欺詐樣本。
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
import pickle
import matplotlib.pyplot as plt
import tensorflow as tf
import seaborn as sns
plt.style.use('seaborn')
from sklearn.model_selection import train_test_split
from keras.models import Model, load_model
from keras.layers import Input, Dense
from keras.callbacks import ModelCheckpoint
from keras import regularizers
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_curve, auc, precision_recall_curve
import sys # 匯入sys模組
# 讀取資料 ：信用卡欺詐資料集地址https://www.kaggle.com/mlg-ulb/creditcardfraud
d = pd.read_csv('creditcard.csv')
# 處理缺失值
d.fillna(d.mean(), inplace=True)
# 檢視樣本比例
num_nonfraud = np.sum(d['Class'] == 0)
num_fraud = np.sum(d['Class'] == 1)
plt.bar(['Fraud', 'non-fraud'], [num_fraud, num_nonfraud], color='dodgerblue')
plt.show()
# 刪除時間列，對Amount進行標準化
data = d.drop(['Time'], axis=1)
data['Amount'] = StandardScaler().fit_transform(data[['Amount']])
# 為無監督新穎點檢測方法，只提取負樣本，並且按照8:2切成訓練集和測試集
mask = (data['Class'] == 0)
X_train, X_test = train_test_split(data[mask], test_size=0.2, random_state=0)
X_train = X_train.drop(['Class'], axis=1).values
X_test = X_test.drop(['Class'], axis=1).values
# 提取所有正樣本，作為測試集的一部分
X_fraud = data[~mask].drop(['Class'], axis=1).values
# 構建Autoencoder網路模型
# 隱藏層節點數分別為16，8，8，16
# epoch為5，batch size為32
input_dim = X_train.shape[1]
encoding_dim = 16
num_epoch = 5
batch_size = 32
input_layer = Input(shape=(input_dim, ))
encoder = Dense(encoding_dim, activation="tanh",
activity_regularizer=regularizers.l1(10e-5))(input_layer)
encoder = Dense(int(encoding_dim / 2), activation="relu")(encoder)
decoder = Dense(int(encoding_dim / 2), activation='tanh')(encoder)
decoder = Dense(input_dim, activation='relu')(decoder)
autoencoder = Model(inputs=input_layer, outputs=decoder)
autoencoder.compile(optimizer='adam',
loss='mean_squared_error',
metrics=['mae'])
# 模型儲存為model.keras，並開始訓練模型
checkpointer = ModelCheckpoint(filepath="model.keras",
verbose=0,
save_best_only=True)
history = autoencoder.fit(X_train, X_train,
epochs=num_epoch,
batch_size=batch_size,
shuffle=True,
validation_data=(X_test, X_test),
verbose=1,
callbacks=[checkpointer]).history
# 讀取模型
autoencoder = load_model('model.keras')
# 利用autoencoder重建測試集
pred_test = autoencoder.predict(X_test)
# 重建欺詐樣本
pred_fraud = autoencoder.predict(X_fraud)
# 計算重構MSE和MAE誤差
mse_test = np.mean(np.power(X_test - pred_test, 2), axis=1)
mse_fraud = np.mean(np.power(X_fraud - pred_fraud, 2), axis=1)
mae_test = np.mean(np.abs(X_test - pred_test), axis=1)
mae_fraud = np.mean(np.abs(X_fraud - pred_fraud), axis=1)
mse_df = pd.DataFrame()
mse_df['Class'] = [0] * len(mse_test) + [1] * len(mse_fraud)
mse_df['MSE'] = np.hstack([mse_test, mse_fraud])
mse_df['MAE'] = np.hstack([mae_test, mae_fraud])
mse_df = mse_df.sample(frac=1).reset_index(drop=True)
# 輸出標籤為Fraud的資料索引
fraud_indices = mse_df[mse_df['Class'] == 1].index
print("Label為Fraud的資料索引：", fraud_indices)
# 分別畫出測試集中正樣本和負樣本的還原誤差MAE和MSE
markers = ['o', '^']
colors = ['dodgerblue', 'coral']
labels = ['Non-fraud', 'Fraud']
plt.figure(figsize=(14, 5))
plt.subplot(121)
for flag in [1, 0]:
temp = mse_df[mse_df['Class'] == flag]
plt.scatter(temp.index,
temp['MAE'],
alpha=0.7,
marker=markers[flag],
c=colors[flag],
label=labels[flag])
plt.title('Reconstruction MAE')
plt.ylabel('Reconstruction MAE'); plt.xlabel('Index')
plt.subplot(122)
for flag in [1, 0]:
temp = mse_df[mse_df['Class'] == flag]
plt.scatter(temp.index,
temp['MSE'],
alpha=0.7,
marker=markers[flag],
c=colors[flag],
label=labels[flag])
plt.legend(loc=[1, 0], fontsize=12); plt.title('Reconstruction MSE')
plt.ylabel('Reconstruction MSE'); plt.xlabel('Index')
plt.show()
# 畫出Precision-Recall曲線
plt.figure(figsize=(14, 6))
for i, metric in enumerate(['MAE', 'MSE']):
plt.subplot(1, 2, i+1)
precision, recall, _ = precision_recall_curve(mse_df['Class'], mse_df[metric])
pr_auc = auc(recall, precision)
plt.title('Precision-Recall curve based on %s\nAUC = %0.2f'%(metric, pr_auc))
plt.plot(recall[:-2], precision[:-2], c='coral', lw=4)
plt.xlabel('Recall'); plt.ylabel('Precision')
plt.show()
| 圖5 | 圖6 |
|---|---|
 |
 |
總結
異常檢測作為一項關鍵技術，在當今多樣化的應用場景中扮演著越來越重要的角色。異常檢測已經滲透到多個行業中。在網路安全領域，它能夠有效識別網路中的惡意活動、攻擊行為和潛在的系統漏洞。金融行業同樣受益於異常檢測技術，它可以識別交易中的欺詐行為、洗錢活動以及其他違規操作。生物醫學領域中，異常檢測技術的應用有助於在醫療檢查中發現異常結果，如腫瘤和疾病等，從而提高診斷的準確性。物流行業也利用異常檢測來識別和處理包裹丟失、延誤或損壞等問題，提升服務質量和效率。
隨著技術的發展，未來的異常檢測將呈現出幾個顯著的發展趨勢。首先，深度學習技術，包括自編碼器和生成對抗網路等，將在異常檢測領域發揮更大的作用，提供更為精準的識別能力。其次，隨著多模態資料的興起，異常檢測將面臨影像、文字、音訊等多種型別資料的綜合分析挑戰，這要求演算法具備更高的適應性和靈活性。此外，提高異常檢測的解釋性，增強演算法的可解釋性和可信度，也將成為未來研究的重點。然而，異常檢測在發展的同時，也面臨著一些挑戰。資料不均衡問題尤為突出，異常資料相對於正常資料的稀缺性，可能導致演算法效能的下降。高維資料的存在同樣對演算法的效能構成了挑戰，需要研究者採用特殊的處理方法來最佳化效能。實時效能的需求也對演算法的計算複雜度和延遲提出了更高的要求，尤其是在需要快速響應的場景中。
參考文獻
- 異常檢測
- 資料異常檢測（多元高斯模型）
- 異常檢測-Anomaly detection
- 8種Python異常檢測演算法總結