python異常值(outlier)檢測實戰:KMeans + PCA + IsolationForest + SVM + EllipticEnvelope

本文轉載自查看原文 2019-11-30 22:03 1499 python/ 異常值檢測/ IsolationForest/ KMeans/ 異常值/ outlier/ 異常值處理abnormal data/ SVM/ EllipticEnvelope/ PCA

機器學習_深度學習_入門經典（博主永久免費教學視頻系列）

https://study.163.com/course/courseMain.htm?courseId=1006390023&share=2&shareId=400000000398149

轉載https://blog.csdn.net/weixin_42608414/article/details/89092501

作者:Susan Li ,原文:https://towardsdatascience.com/time-series-of-price-anomaly-detection-13586cd5ff46

略有增刪

異常值檢測（outlier）是一種數據挖掘過程，用於確定數據集中發現的異常值並確定其出現的詳細信息。當前自動異常檢測至關重要，因為大量數據無法手動標記異常值。自動異常檢測具有廣泛的應用，例如信用卡欺詐檢測，系統健康監測，故障檢測以及傳感器網絡中的事件檢測系統等。今天我們就通過使用python來實現異常值的自動檢測系統的實戰開發。我們將會使用以下技術來實現異常值檢測:

數據

我們的數據kaggle你可以在這里下載。Expedia是全球最大的在線旅行社（OTA,類似我們的攜程網），它每天為數百萬旅行購物者提供搜索服務其中包括用戶在Expedia網站上搜索酒店的相關信息,如國家，地區，房型，價格，入住天數，入住時間等信息。

我們想通過這個數據集來檢測其中價格的異常值。由於Expedia提供的數據集非常大,為了能很好的演示我們的異常值檢測功能，我們將從Expedia數據集中過濾出一個子集，該子集只包含用戶查詢的酒店標間(srch_room_count=1)和酒店所在地為美國(visitor_location_country_id=219)的信息。字段的含義如下:

prop_id：酒店Id
datetime: 用戶查詢的時間
price_usd：價格(美元)
srch_booking_window:從查詢日期開始的酒店住宿天數
srch_saturday_night_bool：如果住宿從周四晚上開始，小於等於4個晚上(必須包含周六晚上)則為1，否則為0

我們會看到同一家酒店,不同的住宿天數，是否包含周六晚，都會導致標間(單間)價格的不同,我們將從中找出價格的異常值。

df = pd.read_csv('./data/expedia_train.csv')
 
#過濾Id為104517的酒店
df = df.loc[df['prop_id'] == 104517]
#過濾標間
df = df.loc[df['srch_room_count'] == 1]
#219表示美國
df = df.loc[df['visitor_location_country_id'] == 219]
df = df[['date_time', 'price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
df.head(10)

我們看一下數據集元數據信息:

我們發現變量date_time的類型不是datetime類型，這會使我們在做可視化的時候出現問題，所以我們要將date_tiem的類型設置為datetime型，接下來我們主要目的是發現價格(price_usd)的異常值,所以我們首先看一下價格的分布情況:

#將date_time的類型設置為datetime
df['date_time'] = pd.to_datetime(df['date_time'])
df = df.sort_values('date_time')
print(df['price_usd'].describe())
df['price_usd'].hist()

我們發現價格的均值是112，但是最大值卻是5584. 這個一個極端的最大值。似乎所有價格數據都小於500，只有一個極端最大值5584。為了我們在后面能找到更多不是極端的異常值，我們先刪除這個極端最大值。

df = df.loc[df['price_usd'] < 5584]
print(df['price_usd'].describe())
df['price_usd'].hist()

刪除價格的極端最大值以后,價格分布基本趨於正常(略微右偏)。

時間序列可視化

下面我們根據時間對價格進行可視化。

df.plot(x='date_time', y='price_usd', figsize=(12,6))
plt.xlabel('時間')
plt.ylabel('價格(美元)')

a = df.loc[df['srch_saturday_night_bool'] == 0, 'price_usd']
b = df.loc[df['srch_saturday_night_bool'] == 1, 'price_usd']
plt.figure(figsize=(10, 6))
plt.hist(a, bins = 50, alpha=0.5, label='不含周六晚上')
plt.hist(b, bins = 50, alpha=0.5, label='含周六晚上')
plt.legend(loc='upper right')
plt.xlabel('價格')
plt.ylabel('數量')
plt.show();

從上面的直方圖可知含周六晚上的(srch_saturday_night_bool=1)的價格均值要大於不含周六晚上的(srch_saturday_night_bool=1)價格均值。含周末的房價略高一些，這應該是合理的。

基於聚類的異常檢測

k-means是一種廣泛使用的聚類算法。它創建了k個具有相似特性的數據組。不屬於這些組的數據實例可能會被標記為異常。在我們開始k-means聚類之前，我們使用elbow方法來確定最佳聚類數量。

data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
n_cluster = range(1, 20)
kmeans = [KMeans(n_clusters=i).fit(data) for i in n_cluster]
scores = [kmeans[i].score(data) for i in range(len(kmeans))]
 
fig, ax = plt.subplots(figsize=(10,6))
ax.plot(n_cluster, scores)
plt.xlabel('聚類集群數')
plt.ylabel('分數')
plt.title('Elbow 曲線')
plt.show()

為了找出合理的距離中心數，我們嘗試盡可能多的聚類中心數（從1個到20個聚類中心），然后我們畫出Elbow曲線，通過觀察Elbow曲線,我們發現當我們的聚類中心數量增加到10個以上時Elbow曲線趨於收斂，因此我們大致可以將聚類中心數定為10.

下面我們將K-means算法的n_clusters設置為10,然后我們將數據進行3D可視化。

X = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
X = X.reset_index(drop=True)
km = KMeans(n_clusters=10)
km.fit(X)
km.predict(X)
labels = km.labels_
 
fig = plt.figure(1, figsize=(7,7))
ax = Axes3D(fig, rect=[0, 0, 0.95, 1], elev=48, azim=134)
ax.scatter(X.iloc[:,0], X.iloc[:,1], X.iloc[:,2],
          c=labels.astype(np.float), edgecolor="k")
ax.set_xlabel("price_usd")
ax.set_ylabel("srch_booking_window")
ax.set_zlabel("srch_saturday_night_bool")
plt.title("K Means", fontsize=14);

接下來我們要確定需要保留數據中的哪些主要成分(特征)

data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
X = data.values
 
#標准化處理,均值為0,標准差為1
X_std = StandardScaler().fit_transform(X)
mean_vec = np.mean(X_std, axis=0)
#計算協方差
cov_mat = np.cov(X_std.T)
 
#計算特征值和特征向量
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
 
#每個特征值對應一組特征向量
eig_pairs = [ (np.abs(eig_vals[i]),eig_vecs[:,i]) for i in range(len(eig_vals))]
eig_pairs.sort(key = lambda x: x[0], reverse= True)
 
#特征值求和
tot = sum(eig_vals)
 
#每個主要成分的解釋方差
var_exp = [(i/tot)*100 for i in sorted(eig_vals, reverse=True)] 
#累計的解釋方差
cum_var_exp = np.cumsum(var_exp) 
 
plt.figure(figsize=(10, 5))
plt.bar(range(len(var_exp)), var_exp, alpha=0.3, align='center', label='獨立的解釋方差', color = 'g')
plt.step(range(len(cum_var_exp)), cum_var_exp, where='mid',label='累積解釋方差')
plt.ylabel('解釋方差率')
plt.xlabel('主成分')
plt.legend(loc='best')
plt.show();

我們首先對數據進行標准化處理(StandardScaler)，然后再計算特征變量之間的協方差矩陣,協方差矩陣反應了特征變量之間的相關性，如果兩個特征變量之間的協方差為正則說明它們之間是正相關關系,如果為負則說明它們之間是負相關關系，如果為0則說明特征變量之間是相互獨立的關系,不存在相關關系(有時候我們也會用相關系數矩陣來代替協方差矩陣)。然后在協方差矩陣的基礎上又計算了協方差矩陣的特征值和特征向量，根據特征值計算出每個主成分(特征)的解釋方差,以及累計解釋方差，我們這樣做的目的是為了下一步做主成分分析(PCA)挑選出特征變量中的主成分。我們挑選前2個主成分，因為它們的累計解釋方差為80%。

從上圖可知我們的三個主成分,第一個主成分(特征)解釋了將近50%的方差變化,第二個主成分解釋了近30%的方差變化,那么前2個主成分解釋了近80%的方差。因此接下來我們將使用PCA算法進行降維並將設置參數n_components=2。

data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
X = data.values
 
#標准化處理,均值為0,標准差為1
X_std = StandardScaler().fit_transform(X)
data = pd.DataFrame(X_std)
 
#將特征維度降到2
pca = PCA(n_components=2)
data = pca.fit_transform(data)
# 降維后將2個新特征進行標准化處理
scaler = StandardScaler()
np_scaled = scaler.fit_transform(data)
data = pd.DataFrame(np_scaled)
 
kmeans = [KMeans(n_clusters=i).fit(data) for i in n_cluster]
df['cluster'] = kmeans[9].predict(data)
df.index = data.index
df['principal_feature1'] = data[0]
df['principal_feature2'] = data[1]
df.head()

基於聚類的異常檢測中的假設是，如果我們對數據進行聚類，則正常數據將屬於聚類，而異常將不屬於任何聚類或屬於小聚類。我們使用以下步驟來查找和可視化異常值。

計算每個數據點與其最近的聚類中心之間的距離。最大的距離被認為是異常的。
設定一個異常值的比例outliers_fraction為1%,這樣設置是因為在標准正太分布的情況下（N(0,1)）我們一般認定3個標准差以外的數據為異常值,3個標准差以內的數據包含了數據集中99%以上的數據,所以剩下的1%的數據可以視為異常值。
根據異常值比例outliers_fraction，計算異常值的數量number_of_outliers
設定一個判定異常值的閾值threshold
通過閾值threshold來判定數據是否為異常值
對數據進行可視化(包含正常數據和異常數據)

# 計算每個數據點到其聚類中心的距離
def getDistanceByPoint(data, model):
    distance = pd.Series()
    for i in range(0,len(data)):
        Xa = np.array(data.loc[i])
        Xb = model.cluster_centers_[model.labels_[i]]
        distance.set_value(i, np.linalg.norm(Xa-Xb))
    return distance
 
#設置異常值比例
outliers_fraction = 0.01
 
# 得到每個點到取聚類中心的距離，我們設置了10個聚類中心，kmeans[9]表示有10個聚類中心的模型
distance = getDistanceByPoint(data, kmeans[9])
 
#根據異常值比例outliers_fraction計算異常值的數量
number_of_outliers = int(outliers_fraction*len(distance))
 
#設定異常值的閾值
threshold = distance.nlargest(number_of_outliers).min()
 
#根據閾值來判斷是否為異常值
df['anomaly1'] = (distance >= threshold).astype(int)
 
#數據可視化
fig, ax = plt.subplots(figsize=(10,6))
colors = {0:'blue', 1:'red'}
ax.scatter(df['principal_feature1'], df['principal_feature2'], c=df["anomaly1"].apply(lambda x: colors[x]))
plt.xlabel('principal feature1')
plt.ylabel('principal feature2')
plt.show();

上圖中紅色的點即是被認定的異常值，它們大約占總數據量的1%。

df = df.sort_values('date_time')
df['date_time_int'] = df.date_time.astype(np.int64)
fig, ax = plt.subplots(figsize=(12,6))
 
a = df.loc[df['anomaly1'] == 1, ['date_time_int', 'price_usd']] #anomaly
 
ax.plot(df['date_time_int'], df['price_usd'], color='blue', label='正常值')
ax.scatter(a['date_time_int'],a['price_usd'], color='red', label='異常值')
plt.xlabel('Date Time Integer')
plt.ylabel('價格(美元)')
plt.legend()
plt.show()

從上圖可知,經過PCA和KMeans計算出的異常值，它們的價格大多位於價格區間的最高點和最低點處，這應該是合理的。

孤立森林(IsolationForest)異常檢測

IsolationForest算法它是一種集成算法(類似於隨機森林)主要用於挖掘異常(Anomaly)數據，或者說離群點挖掘，總之是在一大堆數據中，找出與其它數據的規律不太符合的數據。該算法不采樣任何基於聚類或距離的方法，因此他和那些基於距離的的異常值檢測算法有着根本上的不同，孤立森林認定異常值的原則是異常值是少數的和不同的數據。它通常用於網絡安全中的攻擊檢測和流量異常等分析，金融機構則用於挖掘出欺詐行為。

當我們使用IsolationForest算法時需要設置一個異常值比例的參數contamination，該參數的作用類似於之前的outliers_fraction。
使用 fit 方法對孤立森林模型進行訓練
使用 predict 方法去發現數據中的異常值。返回1表示正常值，-1表示異常值。

data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
scaler = StandardScaler()
np_scaled = scaler.fit_transform(data)
data = pd.DataFrame(np_scaled)
# 訓練孤立森林模型
model =  IsolationForest(contamination=outliers_fraction)
model.fit(data)
 
#返回1表示正常值，-1表示異常值
df['anomaly2'] = pd.Series(model.predict(data)) 
 
fig, ax = plt.subplots(figsize=(10,6))
a = df.loc[df['anomaly2'] == -1, ['date_time_int', 'price_usd']] #異常值
ax.plot(df['date_time_int'], df['price_usd'], color='blue', label = '正常值')
ax.scatter(a['date_time_int'],a['price_usd'], color='red', label = '異常值')
plt.legend()
plt.show();

從上圖可知，使用孤立森林預測的異常值，它們的價格大多位於價格區間的最高點或最低點處。

支持向量機(SVM)的異常檢測

SVM通常應用於監督式學習，但OneClassSVM算法可用於將異常檢測這樣的無監督式學習，它學習一個用於異常檢測的決策函數其主要功能將新數據分類為與訓練集相似的正常值或不相似的異常值。

OneClassSVM

OneClassSVM的思想來源於這篇論文,SVM使用大邊距的方法，它用於異常檢測的主要思想是:將數據密度較高的區域分類為正，將數據密度較低的區域分類為負,如下圖所示:

在訓練OneClassSVM模型時，我們需要設置參數nu = outliers_fraction，它是訓練誤差分數的上限和支持向量分數的下限，並且必須在0和1之間。基本上它代表了我們期望的異常值在我們的數據集中的比例。
指定要在算法中使用的核類型：rbf。它使SVM能夠使用非線性函數將超空間投影到更高維度。
gamma是RBF內核類型的參數，並控制各個訓練樣本的影響 - 這會影響模型的“平滑度”。
predict 對數據進行分類，因為我們的模型是單類模型，所以返回+1或-1，-1表示是異常值，1表示是正常值。

data = df[['price_usd', 'srch_booking_window', 'srch_saturday_night_bool']]
scaler = StandardScaler()
np_scaled = scaler.fit_transform(data)
data = pd.DataFrame(np_scaled)
# 訓練 oneclassSVM 模型
model = OneClassSVM(nu=outliers_fraction, kernel="rbf", gamma=0.01)
model.fit(data)
 
df['anomaly3'] = pd.Series(model.predict(data))
fig, ax = plt.subplots(figsize=(10,6))
 
a = df.loc[df['anomaly3'] == -1, ['date_time_int', 'price_usd']] #anomaly
 
ax.plot(df['date_time_int'], df['price_usd'], color='blue', label ='正常值')
ax.scatter(a['date_time_int'],a['price_usd'], color='red', label = '異常值')
plt.legend()
plt.show();

從上圖可知，使用OneClassSVM預測的異常值，它們的價格大多位於價格區間的最高點或最低點處。

基於高斯概分布的異常檢測

高斯分布也稱為正態分布。它可以被用來進行異常值檢測，不過我們首先要假設我們的數據是正態分布的。不過這個假設不能適應於所有數據集。但如果我們做了這種假設那么它將會有一種有效的方法來發現異常值。

Scikit-Learn的EllipticEnvelope模型，它在假設我們的數據是多元高斯分布的基礎上計算出高斯分布的一些關鍵參數過程。過程大致如下:

根據前面定義的類別創建兩個不同的數據集： search_Sat_night和Search_Non_Sat_night。
在每個類別應用EllipticEnvelope（高斯分布）。
我們設置contamination參數，它表示我們數據集中異常值的比例。
使用decision_function來計算給定數據的決策函數。它等於移位的馬氏距離(Mahalanobis distances)。異常值的閾值為0，這確保了與其他異常值檢測算法的兼容性。
使用predict 來預測數據是否為異常值(1 正常值, -1 異常值)

df_class0 = df.loc[df['srch_saturday_night_bool'] == 0, 'price_usd']
df_class1 = df.loc[df['srch_saturday_night_bool'] == 1, 'price_usd']
 
envelope =  EllipticEnvelope(contamination = outliers_fraction) 
X_train = df_class0.values.reshape(-1,1)
envelope.fit(X_train)
df_class0 = pd.DataFrame(df_class0)
df_class0['deviation'] = envelope.decision_function(X_train)
df_class0['anomaly'] = envelope.predict(X_train)
 
envelope =  EllipticEnvelope(contamination = outliers_fraction) 
X_train = df_class1.values.reshape(-1,1)
envelope.fit(X_train)
df_class1 = pd.DataFrame(df_class1)
df_class1['deviation'] = envelope.decision_function(X_train)
df_class1['anomaly'] = envelope.predict(X_train)
 
df_class = pd.concat([df_class0, df_class1])
df['anomaly5'] = df_class['anomaly']
fig, ax = plt.subplots(figsize=(10, 6))
a = df.loc[df['anomaly5'] == -1, ('date_time_int', 'price_usd')] 
ax.plot(df['date_time_int'], df['price_usd'], color='blue')
ax.scatter(a['date_time_int'],a['price_usd'], color='red')
plt.show()

從上圖可知，使用EllipticEnvelope預測的異常值，它們的價格大多位於價格區間的最高點處在最低點處沒有出現異常值。

到目前為止，我們已經用四種不同的方法進行了價格異常檢測。因為我們的異常檢測是無監督學習。在構建模型之后，我們不知道他們的異常檢測效果怎么樣，因為我們沒有辦法可以對他們進行測試。通常異常檢測只有在實際的應用場景中才能測試出它的效果。

參考

Introduction to Anomaly Detection

sklearn.ensemble.IsolationForest

sklearn.svm.OneClassSVM

sklearn.covariance.EllipticEnvelope

Unsupervised Anomaly Detection | kaggle