LOF離羣因子檢測算法及python3實現

LOF離羣因子檢測算法及python3實現

轉載聲明

本文首發於 知乎專欄-數據科學與python小記html

原文連接 LOF離羣因子檢測算法及python3實現 - Suranyi的文章 - 知乎python

1、相關背景

1.1 異常檢測算法

隨着數據挖掘技術的快速發展,人們在關注數據總體趨勢的同時,開始愈來愈關注那些明顯偏離數據總體趨勢的離羣數據點,由於這些數據點每每蘊含着更加劇要的信息,而處理這些離羣數據要依賴於相應的數據挖掘技術。 離羣點挖掘的目的是有效的識別出數據集中的異常數據,而且挖掘出數據集中有意義的潛在信息。出現離羣點的緣由各有不一樣,其中主要有如下幾種狀況:ios

  • 數據來源於異類:如欺詐、入侵、疾病爆發、不一樣尋常的實驗結果等。這類離羣點的產生機制偏離正常軌道,每每是被關注的重點
  • 數據變量固有的變化:天然、週期發生等反映數據集的分佈特徵,如氣候的忽然變化、顧客新型購買模式、基因突變等
  • 數據測最或收集偏差:主要是系統偏差、人爲的數據讀取誤差、測量設備出現故障或「噪音」干擾
  • 隨機或人爲偏差:主要緣由是實驗平臺存在隨機機制,同時人爲在數據錄入等過程當中可能出現的偏差

離羣點檢測具備很是強的實際意義,在相應的應用領域有着普遍前景。其中工程應用領域主要有如下幾個方面:算法

  • 欺詐檢測:信用卡的不正當行爲,如信用卡、社會保障的欺詐行爲或者是銀行卡、電話卡的欺詐使用等
  • 工業檢測:如計算機網絡的非法訪問等
  • 活動監控:經過實時檢測手機活躍度或者是股權市場的可疑交易,從而實現檢測移動手機詐騙行爲等
  • 網絡性能:計算機網絡性能檢測(穩健性分析),檢測網絡堵塞狀況等
  • 天然生態應用領域:生態系統失調、異常天然氣候的發現等
  • 公共服務領域:公共衛生中的異常疾病的爆發、公共安全中的突發事件的發生等

目前,隨着離羣點檢測技術的日漸成熟,在將來的發展中將會應用在更多的行業當中,而且能更好地爲人類的決策提供指導做用。編程

離羣點檢測的一個目標是從看似雜亂無章的大量數據中挖掘有價值的信息,使這些數據更好地爲咱們的平常生活所服務。可是現實生活中的數據每每具備成百上千的維度,而且數據量極大,這無疑給目前現有的離羣點檢測方法帶來大難題。傳統的離羣點檢測方法雖然在各自特定的應用領域裏表現出很好效果,但在高維大數據集中卻再也不適用。所以如何把離羣點檢測方法有效地應用於大數據、高維度數據,是目前離羣點檢測方法的首要目標之一。數組

1.2 爲何是異常點檢測

分類學習算法在訓練時有一個共同的基本假設:不一樣類別的訓練樣例數目至關。若是不一樣類別的訓練樣例數目稍有差異,一般影響不大,但若差異很大,則會對學習過程形成困擾。 ——周志華《機器學習》安全

2018 年 Mathorcup 數學建模競賽 B 題——「品牌手機目標用戶的精準營銷」中就出現這樣的問題。原題中,檢測用戶數 1萬+,購買手機用戶500+。若使用分類算法,那麼分類器極可能會返回這樣的一個算法:「全部用戶都不會購買這款手機」,分類的正確率高達96%,但顯然沒有實際意義,由於它並不能預測出任何的正例。bash

這類問題稱爲「類別不平衡」,指不一樣類別的訓練樣例數目差異很大的狀況(上例中,購買的與沒有購買用戶數量差異大)。處理這類問題每每採用「欠採樣」、「過採樣」進行數據處理,但經過這樣的方法,可能會損失原始數據中的信息。所以,從離羣點的角度出發,將購買行爲視爲「異常」,進行離羣點挖掘。網絡

1.3 離羣點挖掘方法

離羣點分析

1.4 LOF 算法背景

基於密度的離羣點檢測方法的關鍵步驟在於給每一個數據點都分配一個離散度,其主要思想是:針對給定的數據集,對其中的任意一個數據點,若是在其局部鄰域內的點都很密集,那麼認爲此數據點爲正常數據點,而離羣點則是距離正常數據點最近鄰的點都比較遠的數據點。一般有閾值進行界定距離的遠近。在基於密度的離羣點檢測方法中,最具備表明性的方法是局部離羣因子檢測方法 (Local Outlier Factor, LOF)。app

2、算法簡介

在衆多的離羣點檢測方法中,LOF 方法是一種典型的基於密度的高精度離羣點檢測方法。在 LOF 方法中,經過給每一個數據點都分配一個依賴於鄰域密度的離羣因子 LOF,進而判斷該數據點是否爲離羣點。若 LOF \gg 1, 則該數據點爲離羣點;若 LOF 接近於 1,則該數據點爲正常數據點。

2.1 距離度量尺度

設對於沒有相同點的樣本集合 D ,假設共有 n 個檢測樣本,數據維數爲 m ,對於

\forall X_i=(x_{i1},x_{i2},\cdots,x_{im})\in R\qquad i=1,2,\cdots,n\\

針對數據集 D 中的任意兩個數據點 X_i,X_j ,定義以下幾種經常使用距離度量方式

2.1.1 Eucild(歐幾里得)距離:

\text{Euclid}(X_i,X_j)=\sqrt{\sum^n_{k=1} (X_{ij}-X_{jk})^2}\tag{1}

2.1.2 Hamming(漢明)距離:

\text{Hamming}(X_i,X_j)=\sum^n_{k=1}|X_{ik}\cap X_{jk}|\tag{2}

漢明距離使用在數據傳輸差錯控制編碼裏面,用於度量信息不相同的位數。

X_1=(1,2,3,3,2),X_2=(1,2,4,3,1) ,易見 X_1X_2 中有 2 位數字不相同,所以 X_1X_2 的漢明距離爲 2

對於數據處理,一種技巧是先對連續數據進行分組,化爲分類變量(分組變量),對分類變量能夠引入漢明距離進行度量。——沃茲基 · 碩德

沃茲基 · 碩德

2.1.3 Mahalanobis(馬氏)距離:

設樣本集 D 的協方差矩陣爲 \sum{} ,記其逆矩陣爲 \sum\,^{-1}\sum{ } 可逆,對 \sum{}SVD 分解(奇異值分解),獲得:

\sum=UDV^*\tag{3}

\sum{} 不可逆,則使用廣義逆矩陣 \sum\,^{+}代替 \sum\,^{-1},對其求彭羅斯廣義逆,有:

\sum\,^+=UD^+V^*\tag{4}

則兩個數據點 X_i,X_j 的馬氏距離爲:

M_D(X_i,X_j)=\sqrt{(X_i-Y_i)^T\sum\,^{-1}(X_i-Y_i)}\tag{5}

馬氏距離表示數據的協方差距離,利用 Cholesky 變換處理不一樣維度之間的相關性和度量尺度變換的問題,是一種有效計算樣本集之間的類似度的方法。

2.1.4 球面距離

球面距離實際上是在歐式距離基礎上進行轉換獲得的,並非一種獨特的距離度量方式,在地理信息轉換中常用,本文對此進行詳細介紹。

符號 說明 符號 說明 符號 說明
O 球體球心 \varphi_A A點緯度餘角 O_A A點緯圓圓心
R 球體球徑 \varphi_B B點緯度餘角 O_B B點緯圓圓心
A,B 待測距離點 \psi_{A-B} A,B點經度差 R_A A點緯圓半徑
B' A的緯圓與 B的經圓 交點 R_A B點緯圓半徑

圖:求A、B兩點的球面距離

圖:求A、B兩點的球面距離

設 A,B 兩點的球面座標爲 (x_A,y_A),(x_B,y_B) 。若該球體爲地球,則 x,y 分別表明緯度和經度。(注:下文的 \varphi 爲 x 的餘角,便於推導所使用的記號)

如圖所示,鏈接OA、OB、AB,在 \text{Rt}\Delta OO_AA\text{Rt}\Delta OO_BB 中計算:

\begin{split} O_AA=&R_A=R\cdot \sin\varphi_A\\ O_BB=&R_B=R\cdot \sin\varphi_B\\ O_AO_B=&R\cdot\cos\varphi_A+R\cdot\cos\varphi_B \end{split} \tag{6}

因爲 O_AAO_BB 是異面直線, O_AO_B 是它們的公垂線,所成角度經度差爲 \psi_{A-B} ,利用異面直線上兩點距離公式:

AB^2=O_AA^2+O_BB^2+O_AO_B^2-2\cdot O_AA\cdot O_BB\cdot\cos\psi_{A-B}\tag{7}

\Delta OAB 中,由余弦定理:

\begin{split} \cos\angle AOB&=\frac{OA^2+OB^2-AB^2}{2\cdot OA\cdot OB}\\ &=\sin\varphi_A\sin\varphi_B\cos\psi_{A-B}-\cos\varphi_A\cos\varphi_B \end{split}\tag{8}

因爲此處的 \varphi_A,\varphi_B 表明緯度的補角,對其進行轉換:

\begin{split} \cos\angle AOB&=\sin\varphi_A\sin\varphi_B\cos\psi_{A-B}-\cos\varphi_A\cos\varphi_B\\ &=\cos x_A\cos x_B\cos(y_A-y_B)+\sin x_A\sin x_B \end{split}\tag{9}

所以,點 A,B 的球面距離爲:

\begin{split} d(A,B)&=\frac{\pi R}{180}\cdot\arccos\angle AOB\\ &=\frac{\pi R}{180}\cdot\arccos(\cos x_A\cos x_B\cos(y_A-y_B)+\sin x_A\sin x_B) \end{split}\tag{10}

以上度量方式能夠任意選擇

以上度量方式能夠任意選擇

此外還有 Chebyshev(切比雪夫)距離、Minkowski(閔科夫斯基)距離、絕對值距離、Lance & Williams 距離,具體問題具體分析,選擇合適的度量方式。

統一使用 d(A,B) 表示點 A 和點 B 之間的距離。根據定義,易知交換律成立:

d(A,B)=d(B,A)\tag{11}

2.2 第 k 距離

定義 d_k(O) 爲點 O 的第 k 距離, d_k(O)=d(O,P) ,知足以下條件

  • (1) 在集合中至少存在 k 個點 P'\in D\text{\ }\{O\} ,使得 d(O,P')\le d(O,P)
  • (2) 在集合中至多存在 k-1 個點 P'\in D\text{\ }\{O\} ,使得 d(O,P')\lt d(O,P)

簡言之:點 P 是距離 O 最近的第 k 個點

點P是離O最近的第5個點,第5距離內部有4個點,第5距離內共有6個點。點O的第5距離與第6距離相等

2.3 k 距離鄰域

定義 設 N_{k}(O) 爲點 O 的第 k 距離鄰域,知足:

N_{k}(O)=\{P'\in D\text{\ \{O\}}\mid d(O,P')\le d_k(O)\}\tag{12}

此處的鄰域概念與國內高數教材略有不一樣(具體的點,而非區間)。該集合中包含全部到點 O 距離小於點 O 第 k 鄰域距離的點。易知有N_k(O)\ge k,如上圖,點 O 的第 5 距離鄰域爲:

N_5(O)=\{P_1,P_2,P_3,P_4,P_5,P_6\}\tag{13}

2.4 可達距離

定義 點 P 到點 O 的第 k 可達距離定義爲:

d_k(O,P)=\max\{d_k(O),d(O,P)\}\tag{14}

即點 P 到點 O 的第 k 可達距離至少是點 O 的第 k 距離。距離 O 點最近的 k 個點,它們到 O 的可達距離被認爲是至關的,且都等於 d_k(O)

點 P 到點 O 的第5可達距離,在計算樣本點的可達密度時,此部分老是取d(O)的,即與第k距離有關。若求新樣本點在舊樣本域內的離羣密度,式(14)纔會發揮做用。(如本文 3.1.4 的 preditc,此時要求計算predict每一個點在data下的的密度而不考慮自身數據集影響)

2.5 局部可達密度

定義 局部可達密度定義爲:

\rho_k(O)=\frac{1}{\displaystyle\sum_{P\in N_{k}(O)} d_k(O,P)/k}=\frac{k}{\displaystyle\sum_{P\in N_{k}(O)} d_k(O,P)}\tag{15}

表示點 O 的第 k 鄰域內全部點到 O 的平都可達距離,位於第 k 鄰域邊界上的點即便個數大於 1,也仍將該範圍內點的個數計爲 k 。若是 O 和周圍鄰域點是同一簇,那麼可達距離越可能爲較小的 d_k(O) ,致使可達距離之和越小,局部可達密度越大。若是 O 和周圍鄰域點較遠,那麼可達距離可能會取較大值 d(O,P) ,致使可達距離之和越大,局部可達密度越小。

部分資料這裏使用 |N_k(O)| 而不是 k 。筆者查閱大量資料及數據測試後認爲,此處應爲 k ,不然 \rho_k(O) 會由於過多點在內部同一圓環上(如式13中的 P_5,P_6 位於同一圓環上)而致使 |N_k(O)| 是一個很小的數,提示此處的密度低,可能爲離羣值。
此外,本文 2.1 開頭指出「沒有樣本點重合」在這裏也能獲得解釋:若是考慮重合樣本點,可能會形成此處的可達密度爲 \infty 或下文的 LOF_k\frac{\infty}{\infty} 形式,計算上帶來困擾。

2.6 局部離羣因子

LOF_k(O)=\frac{\displaystyle\sum_{P\in N_{k} (O)}\frac{\rho_k(P)}{\rho_k(O)}}{k}\tag{16}

表示點 O 的鄰域 N_{k } (O) 內其餘點的局部可達密度與點 O 的局部可達密度之比的平均數。若是這個比值越接近1,說明 O 的鄰域點密度差很少, O 可能和鄰域同屬一簇;若是這個比值小於1,說明 O 的密度高於其鄰域點密度, O 爲密集點;若是這個比值大於1,說明 O 的密度小於其鄰域點密度,O 多是異常點。

點O的 LOF 值大於1,提示爲可能的異常點

3、LOF離羣因子檢測算法python3 實現(第 4 部分進行改進)

此部分先介紹 sklearn 提供的函數,再介紹逐步編程實現方法。

3.1 導入 sklearn 模塊實現

在 python3 中,sklearn 模塊提供了 LOF 離羣檢測算法

3.1.1 導入模塊

import pandas as pd
import matplotlib.pyplot as plt
複製代碼

3.1.2 核心函數

clf = LocalOutlierFactor(n_neighbors=20, algorithm='auto', contamination=0.1, n_jobs=-1)
複製代碼
  • n_neighbors = 20:即上文說起的 k ,檢測的鄰域點個數超過樣本數則使用全部的樣本進行檢測
  • algorithm = 'auto':使用的求解算法,使用默認值便可
  • contamination = 0.1:範圍爲 (0, 0.5),表示樣本中的異常點比例,默認爲 0.1
  • n_jobs = -1:並行任務數,設置爲-1表示使用全部CPU進行工做
  • p = 2:距離度量函數,默認使用歐式距離。(其餘距離模型使用較少,這裏不做介紹。具體參考官方文檔
clf.fit(data)
複製代碼

無監督學習,只須要傳入訓練數據data,傳入的數據維度至少是 2 維

clf.kneighbors(data)
複製代碼
  • 做用:獲取第 k 距離鄰域內的每個點到中心點的距離,並按從小到大排序

  • 返回 2\times n\times k 數組,[距離,樣本索引]

-clf._decision_function(data)
複製代碼
  • 做用:獲取每個樣本點的 LOF 值,該函數範圍 LOF 值的相反數,須要取反號
  • clf._decision_function 的輸出方式更爲靈活:若使用 clf._predict(data) 函數,則按照原先設置的 contamination 輸出判斷結果(按比例給出判斷結果,異常點返回-1,非異常點返回1)

3.1.3 封裝函數

localoutlierfactor(data, predict, k)

  • 輸入:訓練樣本,測試樣本, k
  • 輸出:每個測試樣本 LOF 值及對應的第 k 距離

plot_lof(result,method)

  • 輸入:LOF 值、閾值
  • 輸出:以索引爲橫座標的 LOF 值分佈圖

lof(data, predict=None, k=5, method=1, plot = False)

  • 輸入:訓練樣本,測試樣本, k 值,離羣閾值,是否繪製 LOF 圖
  • 輸出:離羣點、正常點分類狀況 沒有輸入測試樣本時,默認測試樣本=訓練樣本
def localoutlierfactor(data, predict, k):
    from sklearn.neighbors import LocalOutlierFactor
    clf = LocalOutlierFactor(n_neighbors=k + 1, algorithm='auto', contamination=0.1, n_jobs=-1)
    clf.fit(data)
    # 記錄 k 鄰域距離
    predict['k distances'] = clf.kneighbors(predict)[0].max(axis=1)
    # 記錄 LOF 離羣因子,作相反數處理
    predict['local outlier factor'] = -clf._decision_function(predict.iloc[:, :-1])
    return predict

def plot_lof(result, method):
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
    plt.rcParams['axes.unicode_minus'] = False  # 用來正常顯示負號
    plt.figure(figsize=(8, 4)).add_subplot(111)
    plt.scatter(result[result['local outlier factor'] > method].index,
                result[result['local outlier factor'] > method]['local outlier factor'], c='red', s=50,
                marker='.', alpha=None,
                label='離羣點')
    plt.scatter(result[result['local outlier factor'] <= method].index,
                result[result['local outlier factor'] <= method]['local outlier factor'], c='black', s=50,
                marker='.', alpha=None, label='正常點')
    plt.hlines(method, -2, 2 + max(result.index), linestyles='--')
    plt.xlim(-2, 2 + max(result.index))
    plt.title('LOF局部離羣點檢測', fontsize=13)
    plt.ylabel('局部離羣因子', fontsize=15)
    plt.legend()
    plt.show()

def lof(data, predict=None, k=5, method=1, plot=False):
    import pandas as pd
    # 判斷是否傳入測試數據,若沒有傳入則測試數據賦值爲訓練數據
    try:
        if predict == None:
            predict = data.copy()
    except Exception:
        pass
    predict = pd.DataFrame(predict)
    # 計算 LOF 離羣因子
    predict = localoutlierfactor(data, predict, k)
    if plot == True:
        plot_lof(predict, method)
    # 根據閾值劃分離羣點與正常點
    outliers = predict[predict['local outlier factor'] > method].sort_values(by='local outlier factor')
    inliers = predict[predict['local outlier factor'] <= method].sort_values(by='local outlier factor')
    return outliers, inliers
複製代碼

3.1.4 實例測試

測試數據2017年全國大學生數學建模競賽B題數據

測試數據有倆份文件,進行三次測試:沒有輸入測試樣本、輸入測試樣本、測試樣本與訓練樣本互換

測試 1 沒有輸入測試樣本狀況:任務密度

數據背景:衆包任務價格制定中,地區任務的密度反映任務的密集程度,從而影響任務的訂價,此處不考慮球面距離誤差(即認爲是同一個平面上的點),如今須要使用一個合理的指標刻畫任務的密集程度。

import numpy as np
import pandas as pd

# 根據文件位置自行修改
posi = pd.read_excel(r'已結束項目任務數據.xls')
lon = np.array(posi["任務gps經度"][:])  # 經度
lat = np.array(posi["任務gps 緯度"][:])  # 緯度
A = list(zip(lat, lon))  # 按照緯度-經度匹配

# 獲取任務密度,取第5鄰域,閾值爲2(LOF大於2認爲是離羣值)
outliers1, inliers1 = lof(A, k=5, method = 2)
複製代碼

給定數據中共有835條數據,設置 LOF 閾值爲 2,輸出17個離羣點信息:

\begin{array}{lrrrr}\hline {\text{index}} & {\text{lat}} & {\text{lon}} & {\text{k distances}} & {\text{local outlier factor}} \\\hline 449 & 22.578728 & 114.493610 & 0.219784 & 2.017148 \\ 507 & 22.861397 & 113.853998 & 0.060414 & 2.019448 \\ 384 & 22.649463 & 113.612985 & 0.120596 & 2.020074 \\ 152 & 23.079789 & 113.429745 & 0.049344 & 2.030277 \\ 419 & 22.894453 & 113.409664 & 0.048589 & 2.036357 \\ 405 & 22.842245 & 113.345077 & 0.070277 & 2.052435 \\ 643 & 22.969205 & 114.174155 & 0.114833 & 2.065984 \\ 410 & 22.912291 & 113.350764 & 0.032169 & 2.088684 \\ 409 & 22.839969 & 113.349897 & 0.072327 & 2.102734 \\ 294 & 23.563411 & 113.580385 & 0.018121 & 2.507724 \\ 744 & 23.221139 & 112.924846 & 0.058768 & 2.689194 \\ 224 & 23.338871 & 113.111065 & 0.108733 & 2.869059 \\ 425 & 23.454574 & 113.498304 & 0.132545 & 7.370212 \\ 296 & 23.878398 & 113.539711 & 0.324031 & 15.953839 \\ 295 & 23.625476 & 113.431396 & 0.178774 & 20.452951 \\ 297 & 23.723118 & 113.739427 & 0.222326 & 28.793106 \\ 302 & 23.816108 & 113.957929 & 0.443798 & 29.319828 \\\hline \end{array} \\

繪製數據點檢測狀況分佈以下圖所示,其中藍色表示任務分佈狀況,紅色範圍表示 LOF 值大小:

k=5時檢測狀況:紅色圈越大,LOF 值越大。從圖中能夠看出檢測效果顯著

調整 k 值進行屢次檢測效果,k值越大精度越高

# 繪圖程序
import matplotlib.pyplot as plt
for k in [3,5,10]:
    plt.figure('k=%d'%k)
    outliers1, inliers1 = lof(A, k=k, method = 2)
    plt.scatter(np.array(A)[:,0],np.array(A)[:,1],s = 10,c='b',alpha = 0.5)
    plt.scatter(outliers1[0],outliers1[1],s = 10+outliers1['local outlier factor']*100,c='r',alpha = 0.2)
    plt.title('k=%d' % k)
複製代碼

測試 2 有輸入測試樣本狀況:任務對會員的密度

數據背景:衆包任務價格制定中,地區任務的密度反映任務的密集程度、會員密度反映會員的密集程度。而任務對會員的密度則能夠用於刻畫任務點周圍會員的密集程度,從而體現任務被完成的相對機率。此時訓練樣本爲會員密度,測試樣本爲任務密度。

import numpy as np
import pandas as pd

posi = pd.read_excel(r'已結束項目任務數據.xls')
lon = np.array(posi["任務gps經度"][:])  # 經度
lat = np.array(posi["任務gps 緯度"][:])  # 緯度
A = list(zip(lat, lon))  # 按照緯度-經度匹配

posi = pd.read_excel(r'會員信息數據.xlsx')
lon = np.array(posi["會員位置(GPS)經度"][:])  # 經度
lat = np.array(posi["會員位置(GPS)緯度"][:])  # 緯度
B = list(zip(lat, lon))  # 按照緯度-經度匹配

# 獲取任務對會員密度,取第5鄰域,閾值爲2(LOF大於2認爲是離羣值)
outliers2, inliers2 = lof(B, A, k=5, method=2)
複製代碼

給定訓練樣本中共有1877條數據,測試樣本中共有835條數據,設置 LOF 閾值爲 2,輸出34個離羣點信息:

\begin{array}{lrrrr}\hline {\text{index}} & {\text{lat}} & {\text{lon}} & {\text{k distances}} & {\text{local outlier factor}} \\ \hline 511 & 23.000649 & 113.615440 & 0.065718 & 2.000004 \\ 530 & 22.995226 & 113.800023 & 0.041281 & 2.000528 \\ 89 & 22.704665 & 113.972690 & 0.027174 & 2.025587 \\ 439 & 22.526199 & 113.989221 & 0.026114 & 2.028834 \\ 425 & 23.454574 & 113.498304 & 0.126509 & 2.042262 \\ 55 & 22.776408 & 114.309588 & 0.047912 & 2.042737 \\ 458 & 22.805179 & 113.562852 & 0.110320 & 2.058531 \\ 710 & 22.841581 & 113.983753 & 0.072645 & 2.070924 \\ 215 & 23.374776 & 113.178885 & 0.035750 & 2.091507 \\ 295 & 23.625476 & 113.431396 & 0.171411 & 2.152127 \\ 224 & 23.338871 & 113.111065 & 0.105294 & 2.182644 \\ 152 & 23.079789 & 113.429745 & 0.035729 & 2.183811 \\ 198 & 22.774034 & 113.563485 & 0.112543 & 2.187064 \\ 208 & 22.751494 & 113.583011 & 0.115377 & 2.194350 \\ 206 & 22.751557 & 113.582301 & 0.115689 & 2.198175 \\ 177 & 22.750343 & 113.583523 & 0.116110 & 2.199548 \\ 199 & 22.760969 & 113.567052 & 0.116878 & 2.235774 \\ 633 & 22.998891 & 113.620643 & 0.061624 & 2.242778 \\ 586 & 22.824523 & 112.683258 & 0.213114 & 2.245389 \\ 241 & 23.406491 & 113.169638 & 0.046772 & 2.246444 \\ 106 & 22.670522 & 113.923561 & 0.019000 & 2.259844 \\ 17 & 22.498190 & 113.898482 & 0.035912 & 2.409831 \\ 8 & 22.500012 & 113.895661 & 0.037886 & 2.469827 \\ 403 & 23.013350 & 113.395323 & 0.035931 & 2.934989 \\ 362 & 23.012180 & 113.397405 & 0.037304 & 3.037904 \\ 258 & 23.014037 & 113.399595 & 0.035749 & 3.054759 \\ 128 & 23.168789 & 113.419960 & 0.029593 & 3.078636 \\ 193 & 23.162819 & 113.418328 & 0.028216 & 3.169111 \\ 449 & 22.578728 & 114.493610 & 0.174552 & 3.325349 \\ 221 & 22.710230 & 113.552891 & 0.164642 & 3.344975 \\ 297 & 23.723118 & 113.739427 & 0.223941 & 4.221785 \\ 296 & 23.878398 & 113.539711 & 0.327200 & 4.431575 \\ 384 & 22.649463 & 113.612985 & 0.178139 & 4.708909 \\ 302 & 23.816108 & 113.957929 & 0.445676 & 8.867589 \\\hline \end{array}

繪製數據點檢測狀況分佈以下圖所示,其中藍色表示任務分佈狀況,綠色表示會員分佈狀況,紅色範圍表示 LOF 值大小。

k=5時檢測狀況:紅色圈越大,LOF 值越大。從圖中能夠看出檢測效果顯著,但誤判也較爲嚴重

k=1時檢測狀況,檢測狀況不佳。由於k=1時只對點間的距離進行比較,沒有太大實際意義

# 繪圖程序
import matplotlib.pyplot as plt
for k,v in ([1,5],[5,2]):
    plt.figure('k=%d'%k)
    outliers2, inliers2 = lof(B, A, k=k, method=v)
    plt.scatter(np.array(A)[:,0],np.array(A)[:,1],s = 10,c='b',alpha = 0.5)
    plt.scatter(np.array(B)[:,0],np.array(B)[:,1],s = 10,c='green',alpha = 0.3)
    plt.scatter(outliers2[0],outliers2[1],s = 10+outliers2['local outlier factor']*100,c='r',alpha = 0.2)
    plt.title('k = %d, method = %g' % (k,v))
複製代碼

測試 3 測試樣本與訓練樣本互換:會員對任務的密度

數據背景:衆包任務價格制定中,地區任務的密度反映任務的密集程度、會員密度反映會員的密集程度。而任務對會員的密度則能夠用於刻畫會員周圍任務的密集程度,從而體現會員能接到任務的相對機率。此時訓練樣本爲任務密度,測試樣本爲會員密度。

import numpy as np
import pandas as pd

posi = pd.read_excel(r'已結束項目任務數據.xls')
lon = np.array(posi["任務gps經度"][:])  # 經度
lat = np.array(posi["任務gps 緯度"][:])  # 緯度
A = list(zip(lat, lon))  # 按照緯度-經度匹配

posi = pd.read_excel(r'會員信息數據.xlsx')
lon = np.array(posi["會員位置(GPS)經度"][:])  # 經度
lat = np.array(posi["會員位置(GPS)緯度"][:])  # 緯度
B = list(zip(lat, lon))  # 按照緯度-經度匹配

# 獲取會員對任務密度,取第5鄰域,閾值爲5(LOF大於5認爲是離羣值)
outliers3, inliers3 = lof(A, B, k=5, method=5)
複製代碼

給定訓練樣本中共有835條數據,測試樣本中共有1877條數據,設置 LOF 閾值爲 5,輸出20個離羣點信息:

\begin{array}{lrrrr}\hline {\text{index}} & {\text{lat}} & {\text{lon}} & {\text{k distances}} & {\text{local outlier factor}} \\ \hline 32&22.993463&114.728544 & 0.477421& 5.136386\\ 1779 & 23.587850 & 113.619693& 0.042103 & 6.361438\\ 1385 & 23.563300 & 113.475857 & 0.118328 & 10.011952 \\ 1821 & 22.800504 & 115.374799 & 1.055900 & 10.017484 \\ 1263 & 23.585306 & 113.524752 & 0.077290 & 10.754401 \\ 732 & 23.508004 & 113.686901 & 0.100989 & 13.384237 \\ 1726 & 21.533320 & 111.229119 & 2.162833 & 16.053346 \\ 471 & 21.498823 & 111.106315 & 2.280257 & 16.939837 \\ 38 & 21.679227 & 110.922443 & 2.327916 & 17.287809 \\ 790 & 23.639637 & 113.685956 & 0.124816 & 18.126994 \\ 139 & 23.626787 & 113.436410 & 0.174738 & 20.057431 \\ 81 & 21.202247 & 110.417157 & 3.011823 & 22.475269 \\ 905 & 23.414427 & 113.717294 & 0.183829 & 23.964102 \\ 1707 & 24.293405 & 116.122526 & 2.367451 & 25.220603 \\ 47 & 20.335061 & 110.178827 & 3.743817 & 28.044233 \\ 135 & 24.804130 & 113.605786 & 1.243749 & 36.022419 \\ 21 & 27.124487 & 111.017906 & 4.238630 & 308.053190 \\ 4 & 33.652050 & 116.970470 & 10.642993 & 336.594272 \\ 6 & 29.560903 & 106.239083 & 9.220659 & 539.377888 \\ 1174 & 113.131483 & 23.031824 & 127.133142 & 9272.946682 \\ \hline \end{array}

繪製數據點檢測狀況分佈以下圖所示,其中藍色表示會員分佈狀況,綠色表示任務分佈狀況,紅色範圍表示 LOF 值大小。

k=5時檢測狀況:紅色圈越大,LOF 值越大。從圖中能夠看出檢測效果顯著

# 繪圖程序
plt.figure('k=5')
outliers3, inliers3 = lof(A, B, k=5, method=5)
plt.scatter(np.array(B)[:, 0], np.array(B)[:, 1], s=10, c='b', alpha=0.5)
plt.scatter(np.array(A)[:, 0], np.array(A)[:, 1], s=10, c='green', alpha=0.3)
plt.scatter(outliers3[0], outliers3[1], s=10 + outliers3['local outlier factor'] * 20, c='r', alpha=0.2)
plt.title('k = 5, method = 5')
複製代碼

將 method 設置爲 0 就能輸出每個點的 LOF 值,做爲密度指標。

3.2 逐步編程實現方法

3.2.1 距離度量尺度

distances(A, B,model = 'euclidean')

  • 調用 geopy 模塊中的函數進行球面距離計算
  • 輸入:A,B兩點座標,距離模式('euclidean', 'geo' 分別爲歐式距離和球面距離)
  • 輸出:若A,B爲兩個座標點,則返回座標點距離;若A,B爲矩陣,計算規則以下:
\text{Input:}A=\left[\begin{array}{cc} x_{A1}&y_{A1}\\ x_{A2}&y_{A2}\\ \vdots&\vdots\\ x_{An}&y_{An}\\ \end{array}\right]=\left[\begin{array}{c} P_{A1}\\ P_{A2}\\ \vdots\\ P_{An}\\ \end{array}\right] , B=\left[\begin{array}{cc} x_{B1}&y_{B2}\\ x_{B2}&y_{B2}\\ \vdots&\vdots\\ x_{Bm}&y_{Bm}\\ \end{array}\right]=\left[\begin{array}{c} P_{B1}\\ P_{B2}\\ \vdots\\ P_{Bm}\\ \end{array}\right]\\\tag{17}
\text{Output:}d(A,B)=\left[\begin{array}{cc} d(P_{A1},P_{B1})&d(P_{A1},P_{B2})&\cdots&d(P_{A1},P_{Bm})\\ d(P_{A2},P_{B1})&d(P_{A2},P_{B2})&\cdots&d(P_{A2},P_{Bm})\\ \vdots&\vdots&\ddots&\vdots\\ d(P_{An},P_{B1})&d(P_{An},P_{B2})&\cdots&d(P_{An},P_{Bm})\\ \end{array}\right]\tag{18}
def distances(A, B,model = 'euclidean'):
    '''LOF中定義的距離,默認爲歐式距離,也提供球面距離'''
    import numpy as np
    A = np.array(A); B = np.array(B)
    if model == 'euclidean':
        from scipy.spatial.distance import pdist, squareform
        distance = squareform(pdist(np.vstack([A, B])))[:A.shape[0],A.shape[0]:]
    if model == 'geo':
        from geopy.distance import great_circle
        distance = np.zeros(A.shape[0]*B.shape[0]).reshape(A.shape[0],B.shape[0])
        for i in range(len(A)):
            for j in range(len(B)):
                distance[i,j] = great_circle(A[i], B[j]).kilometers
    if distance.shape == (1,1):
        return distance[0][0]
    return distance
複製代碼

3.2.2 k 鄰域半徑及鄰域點

k_distance(k, instance_A, instance_B, result, model)

  • 調用 3.2.1 的 distances 函數
  • 輸入: k 值,A,B兩點座標,pandas.DataFrame 表用於存儲中間信息量,距離模式
  • 輸出:pandas.DataFrame 表(距離、鄰域點座標)
def k_distance(k, instance_A, instance_B, result, model):
    '''計算k距離鄰域半徑及鄰域點'''
    distance_all = distances(instance_B, instance_A, model)
    # 對 instance_A 中的每個點進行遍歷
    for i,a in enumerate(instance_A):
        distances = {}
        distance = distance_all[:,i]
        # 記錄 instance_B 到 instance_A 每個點的距離,不重複記錄
        for j in range(distance.shape[0]):
            if distance[j] in distances.keys():
                if instance_B[j].tolist() in distances[distance[j]]:
                    pass
                else:
                    distances[distance[j]].append(instance_B[j].tolist())
            else:
                distances[distance[j]] = [instance_B[j].tolist()]
        # 距離排序
        distances = sorted(distances.items())
        if distances[0][0] == 0:
            distances.remove(distances[0])
        neighbours = [];k_sero = 0;k_dist = None
        # 截取前 k 個點
        for dist in distances:
            k_sero += len(dist[1])
            neighbours.extend(dist[1])
            k_dist = dist[0]
            if k_sero >= k:
                break
        # 輸出信息
        result.loc[str(a.tolist()),'k_dist'] = k_dist
        result.loc[str(a.tolist()),'neighbours'] = str(neighbours)
    return result
複製代碼

3.2.3 局部可達密度

local_reachability_density(k,instance_A,instance_B,result, model)

  • 調用 3.2.2 的 k_distance 函數
  • 輸入: k 值,A,B兩點座標,pandas.DataFrame 表用於存儲中間信息量,距離模式
  • 輸出:pandas.DataFrame 表(距離、鄰域點座標、局部可達密度)
def local_reachability_density(k,instance_A,instance_B,result, model):
    '''局部可達密度'''
    result = k_distance(k, instance_A, instance_B, result, model)
    # 對 instance_A 中的每個點進行遍歷
    for a in instance_A:
        # 獲取k_distance中獲得的鄰域點座標,解析爲點座標(字符串 -> 數組 -> 點)
        try:
            (k_distance_value, neighbours) = result.loc[str(a.tolist())]['k_dist'].mean(),eval(result.loc[str(a.tolist())]['neighbours'])
        except Exception:
            (k_distance_value, neighbours) = result.loc[str(a.tolist())]['k_dist'].mean(), eval(result.loc[str(a.tolist())]['neighbours'].values[0])
        # 計算局部可達距離
        reachability_distances_array = [0]*len(neighbours)
        for j, neighbour in enumerate(neighbours):
            reachability_distances_array[j] = max([k_distance_value, distances([a], [neighbour],model)])
        sum_reach_dist = sum(reachability_distances_array)
        # 計算局部可達密度,並儲存結果
        result.loc[str(a.tolist()),'local_reachability_density'] = k / sum_reach_dist
    return result
複製代碼

3.2.4 局部離羣因子

k_distance(k, instance_A, instance_B, result, model)

  • 調用 3.2.3 的 local_reachability_density 函數
  • 輸入: k 值,A,B兩點座標,pandas.DataFrame 表用於存儲中間信息量,距離模式
  • 輸出:pandas.DataFrame 表(距離、鄰域點座標、局部可達密度、離羣因子)
def local_outlier_factor(k,instance_A,instance_B,model):
    '''局部離羣因子'''
    result = local_reachability_density(k,instance_A,instance_B,pd.DataFrame(index=[str(i.tolist()) for i in instance_A]), model)
    # 判斷:若測試數據=樣本數據
    if np.all(instance_A == instance_B):
        result_B = result
    else:
        result_B = local_reachability_density(k, instance_B, instance_B, k_distance(k, instance_B, instance_B, pd.DataFrame(index=[str(i.tolist()) for i in instance_B]), model), model)
    for a in instance_A:
        try:
            (k_distance_value, neighbours, instance_lrd) = result.loc[str(a.tolist())]['k_dist'].mean(),np.array(eval(result.loc[str(a.tolist())]['neighbours'])),result.loc[str(a.tolist())]['local_reachability_density'].mean()
        except Exception:
            (k_distance_value, neighbours, instance_lrd) = result.loc[str(a.tolist())]['k_dist'].mean(), np.array(eval(result.loc[str(a.tolist())]['neighbours'].values[0])), result.loc[str(a.tolist())]['local_reachability_density'].mean()
        finally:
            lrd_ratios_array = [0]* len(neighbours)
            for j,neighbour in enumerate(neighbours):
                neighbour_lrd = result_B.loc[str(neighbour.tolist())]['local_reachability_density'].mean()
                lrd_ratios_array[j] = neighbour_lrd / instance_lrd
            result.loc[str(a.tolist()), 'local_outlier_factor'] = sum(lrd_ratios_array) / k
    return result
複製代碼

3.2.5 封裝函數

lof(k,instance_A,instance_B,k_means=False,$n_clusters$=False,k_means_pass=3,method=1,model = 'euclidean'

  • 調用 3.2.4 的 k_distance 函數
  • 輸入: k 值,A,B兩點座標,是否使用聚類算法,聚類簇數,跳過聚類樣本數低於3的簇,閾值,距離模式
  • 輸出:離羣點信息、正常點信息、LOF 分佈圖
# 函數中的 k_means將在第4部分介紹
def lof(k, instance_A, instance_B, k_means=False, $n_clusters$=False, k_means_pass=3, method=1, model='euclidean'):
    '''A做爲定點,B做爲動點'''
    import numpy as np
    instance_A = np.array(instance_A);
    instance_B = np.array(instance_B)
    if np.all(instance_A == instance_B):
        if k_means == True:
            if $n_clusters$ == True:
                $n_clusters$ = elbow_method(instance_A, maxtest=10)
            instance_A = kmeans(instance_A, $n_clusters$, k_means_pass)
            instance_B = instance_A.copy()
    result = local_outlier_factor(k, instance_A, instance_B, model)
    outliers = result[result['local_outlier_factor'] > method].sort_values(by='local_outlier_factor', ascending=False)
    inliers = result[result['local_outlier_factor'] <= method].sort_values(by='local_outlier_factor', ascending=True)
    plot_lof(result, method) # 該函數見3.1.3
    return outliers, inliers
複製代碼

4、LOF 算法相關思考及改進

4.1 一維數據可使用 LOF 嗎?

本文 3.1 中 sklearn 模塊提供的 LOF 方法進行訓練時會進行數據類型判斷,若數據類型爲list、tuple、numpy.array 則要求傳入數據的維度至少是 2 維。實際上要篩選 1 維數據中的離羣點,直接在座標系中繪製出圖像進行閾值選取判斷也很方便。但此情形下若要使用 LOF 算法,能夠爲數據添加虛擬維度,並賦相同的值:

# data 是原先的1維數據,經過下面的方法轉換爲2維數據
data = list(zip(data, np.zeros_like(data)))
複製代碼

此外,也能夠經過將數據轉化爲 pandas.DataFrame 形式避免上述問題:

data = pd.DataFrame(data)
複製代碼

4.2 多大的 LOF 纔是離羣值?

LOF 計算結果對於多大的值定義爲離羣值沒有明確的規定。在一個平穩數據集中,可能 1.1 已是一個異常值,而在另外一個具備強烈數據波動的數據集中,即便 LOF 值爲 2 可能還是一個正常值。因爲方法的侷限性,數據集中的異常值界定可能存在差別,筆者認爲可使用統計分佈方法做爲參考,再結合數據狀況最終肯定閾值。

基於統計分佈的閾值劃分

將 LOF 異常值分數歸一化到 [0, 1] 區間,運用統計方法進行劃分下面提供使用箱型圖進行界定的方法,根據異常輸出狀況參考選取。

box(data, legend=True)

  • 輸入:LOF 異常分數值,在箱型圖中繪製異常數據(設置爲False不繪製)
  • 輸出:異常識別狀況
def box(data, legend=True):
    import matplotlib.pyplot as plt
    import pandas as pd
    plt.rcParams['axes.unicode_minus'] = False
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.style.use("ggplot")
    plt.figure()
    # 若是不是DataFrame格式,先進行轉化
    if type(data) != pd.core.frame.DataFrame:
        data = pd.DataFrame(data)
    p = data.boxplot(return_type='dict')
    warming = pd.DataFrame()
    y = p['fliers'][0].get_ydata()
    y.sort()
    for i in range(len(y)):
        if legend == True:
            plt.text(1, y[i] - 1, y[i], fontsize=10, color='black', ha='right')
        if y[i] < data.mean()[0]:
            form = '低'
        else:
            form = '高'
        warming = warming.append(pd.Series([y[i], '偏' + form]).T, ignore_index=True)
    print(warming)
    plt.show()
複製代碼

box 函數能夠插入封裝函數 lof 中,傳入 data = predict['local outlier factor'] 實現;也能夠先隨機指定一個初始閾值,(輸出的離羣點、正常點分別命名爲outliers, inliers)再輸入:

box(outliers['local outlier factor'].tolist()+inliers['local outlier factor'].tolist(), legend=True)
複製代碼

此時,交互控制檯中輸出狀況以下左圖所示,箱型圖以下右圖所示。輸出狀況提示咱們從數據分佈的角度上,能夠將 1.4 做爲離羣識別閾值,但實際上取 7 更爲合適(從 2 到 7 間有明顯的斷層,而上文中設定爲 5 是通過屢次試驗後選取的數值)。

設置 legend=False 能夠關閉右圖的標籤

4.3 數據維度過大,還能使用 LOF 算法嗎?

數據維度過大一方面會增大量綱的影響,另外一方面增大計算難度。此時直接使用距離度量的表達形式不合理,並有人爲放大較爲分散數據影響的風險。一種處理方式是採用馬氏距離做爲距離度量的方式(去量綱化)。另外一種處理方式,參考隨機森林的決策思想,能夠考慮在多個維度投影上使用 LOF 並將結果結合起來,以提升高維數據的檢測質量。

集成學習

集成學習:經過構建並結合多個學習器來完成學習任務,一般能夠得到比單一學習器更顯著優越的泛化性能。 ——周志華《機器學習》

數據檢測中進行使用的數據應該是有意義的數據,這就須要進行簡單的特徵篩選,不然不管多麼「離羣」的樣本,可能也沒有多大的實際意義。根據集成學習的思想,須要將數據按維度拆分,對於同類型的數據,這裏假設你已經作好了規約處理(如位置座標能夠放在一塊兒做爲一個特徵「距離」進行考慮),而且數據的維度大於 1,不然使用 4.1 中的數據變換及通常形式 LOF 便可處理。

4.3.1 投票表決模式

投票表決模式認爲每個維度的數據都是同等重要,單獨爲每一個維度數據設置 LOF 閾值並進行比對,樣本的 LOF 值超過閾值則異常票數積 1 分,最終超過票數閾值的樣本認爲是離羣樣本。

localoutlierfactor(data, predict, k)

  • 輸入:訓練樣本,測試樣本, k 值
  • 輸出:每個測試樣本 LOF 值及對應的第 k 距離

plot_lof(result,method)

  • 輸入:LOF 值、閾值
  • 輸出:以索引爲橫座標的 LOF 值分佈圖

ensemble_lof(data, predict=None, k=5, groups=[], method=1, vote_method = 'auto')

  • 輸入:訓練樣本,測試樣本, k 值,組合特徵索引,離羣閾值,票數閾值 組合特徵索引(列位置 - 1):如第一列數據與第二列數據做爲同類型數據,則傳入groups = [[0, 1]]
    • 離羣閾值:每個特徵的離羣閾值,缺乏位數則以 1 替代
    • 票數閾值:正常點離羣因子得票數上限
  • 輸出:離羣點、正常點分類狀況 沒有輸入測試樣本時,默認測試樣本=訓練樣本
def localoutlierfactor(data, predict, k, group_str):
    from sklearn.neighbors import LocalOutlierFactor
    clf = LocalOutlierFactor(n_neighbors=k + 1, algorithm='auto', contamination=0.1, n_jobs=-1)
    clf.fit(data)
    # 記錄 LOF 離羣因子,作相反數處理
    predict['local outlier factor %s' % group_str] = -clf._decision_function(predict.iloc[:, eval(group_str)])
    return predict

def plot_lof(result, method, group_str):
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
    plt.rcParams['axes.unicode_minus'] = False  # 用來正常顯示負號
    plt.figure('local outlier factor %s' % group_str)
    try:
        plt.scatter(result[result > method].index,
                    result[result > method], c='red', s=50,
                    marker='.', alpha=None,
                    label='離羣點')
    except Exception:
        pass
    try:
        plt.scatter(result[result <= method].index,
                    result[result <= method], c='black', s=50,
                    marker='.', alpha=None, label='正常點')
    except Exception:
        pass
    plt.hlines(method, -2, 2 + max(result.index), linestyles='--')
    plt.xlim(-2, 2 + max(result.index))
    plt.title('LOF局部離羣點檢測', fontsize=13)
    plt.ylabel('局部離羣因子', fontsize=15)
    plt.legend()
    plt.show()

def ensemble_lof(data, predict=None, k=5, groups=[], method=1, vote_method = 'auto'):
    import pandas as pd
    import numpy as np
    # 判斷是否傳入測試數據,若沒有傳入則測試數據賦值爲訓練數據
    try:
        if predict == None:
            predict = data.copy()
    except Exception:
        pass
    data = pd.DataFrame(data); predict = pd.DataFrame(predict)
    # 數據標籤分組,默認獨立自成一組
    for i in range(data.shape[1]):
        if i not in pd.DataFrame(groups).values:
            groups += [[i]]
    # 擴充閾值列表
    if type(method) != list:
        method = [method]
        method += [1] * (len(groups) - 1)
    else:
        method += [1] * (len(groups) - len(method))
    vote = np.zeros(len(predict))
    # 計算LOF離羣因子並根據閾值進行票數統計
    for i in range(len(groups)):
        predict = localoutlierfactor(pd.DataFrame(data).iloc[:, groups[i]], predict, k, str(groups[i]))
        plot_lof(predict.iloc[:, -1], method[i], str(groups[i]))
        vote += predict.iloc[:, -1] > method[i]
    # 根據票數閾值劃分離羣點與正常點
    predict['vote'] = vote
    if vote_method == 'auto':
        vote_method = len(groups)/2
    outliers = predict[vote > vote_method].sort_values(by='vote')
    inliers = predict[vote <= vote_method].sort_values(by='vote')
    return outliers, inliers
複製代碼

測試 4 仍然使用測試3的狀況進行分析,此時將經度、緯度設置爲獨立的特徵,分別對兩個維度數據進行識別(儘管單獨的緯度、經度數據沒有太大的實際意義)

import numpy as np
import pandas as pd

posi = pd.read_excel(r'已結束項目任務數據.xls')
lon = np.array(posi["任務gps經度"][:])  # 經度
lat = np.array(posi["任務gps 緯度"][:])  # 緯度
A = list(zip(lat, lon))  # 按照緯度-經度匹配

posi = pd.read_excel(r'會員信息數據.xlsx')
lon = np.array(posi["會員位置(GPS)經度"][:])  # 經度
lat = np.array(posi["會員位置(GPS)緯度"][:])  # 緯度
B = list(zip(lat, lon))  # 按照緯度-經度匹配

# 獲取會員對任務密度,取第5鄰域,閾值分別爲 1.5,2,得票數超過 1 的認爲是異常點 
outliers4, inliers4 = ensemble_lof(A, B, k=5, method=[1.5,2], vote_method = 1)

# 繪圖程序
plt.figure('投票集成 LOF 模式')
plt.scatter(np.array(B)[:, 0], np.array(B)[:, 1], s=10, c='b', alpha=0.5)
plt.scatter(np.array(A)[:, 0], np.array(A)[:, 1], s=10, c='green', alpha=0.3)
plt.scatter(outliers4[0], outliers4[1], s=10 + 1000, c='r', alpha=0.2)
plt.title('k = 5, method = [1.5, 2]')
複製代碼

仍能較爲準確地識別異常點

4.3.2 LOF 異常分數加權模式

異常分數加權模式則是對各維度數據的 LOF 值進行加權,獲取最終的 LOF 得分做爲總體數據的 LOF 得分。權重能夠認爲是特徵的重要程度,也能夠認爲是數據分佈的相對離散程度,若視爲後面一種情形,能夠根據熵權法進行設定,關於熵權法的介紹詳見筆者另外一篇博文

LOF_i=\sum^m_{j=1} w_iLOF_{ij}\tag{19}\\

式中 LOF_{ij} 表示第 i 個數據的第 j 維度 LOF 異常分數值。

localoutlierfactor(data, predict, k)

  • 輸入:訓練樣本,測試樣本, k 值
  • 輸出:每個測試樣本 LOF 值及對應的第 k 距離

plot_lof(result,method)

  • 輸入:LOF 值、閾值
  • 輸出:以索引爲橫座標的 LOF 值分佈圖

ensemble_lof(data, predict=None, k=5, groups=[], method=2, weight=1)

  • 輸入:訓練樣本,測試樣本, k 值,組合特徵索引,離羣閾值,特徵權重 組合特徵索引(列位置 - 1):如第一列數據與第二列數據做爲同類型數據,則傳入groups = [[0, 1]]
    • 離羣閾值:加權 LOF 離羣閾值,默認爲 2
    • 特徵權重:爲每一個維度的 LOF 設置權重,缺乏位數則以 1 代替
  • 輸出:離羣點、正常點分類狀況
  • 沒有輸入測試樣本時,默認測試樣本=訓練樣本
def localoutlierfactor(data, predict, k, group_str):
    from sklearn.neighbors import LocalOutlierFactor
    clf = LocalOutlierFactor(n_neighbors=k + 1, algorithm='auto', contamination=0.1, n_jobs=-1)
    clf.fit(data)
    # 記錄 LOF 離羣因子,作相反數處理
    predict['local outlier factor %s' % group_str] = -clf._decision_function(predict.iloc[:, eval(group_str)])
    return predict

def plot_lof(result, method):
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
    plt.rcParams['axes.unicode_minus'] = False  # 用來正常顯示負號
    plt.scatter(result[result > method].index,
                result[result > method], c='red', s=50,
                marker='.', alpha=None,
                label='離羣點')
    plt.scatter(result[result <= method].index,
                result[result <= method], c='black', s=50,
                marker='.', alpha=None, label='正常點')
    plt.hlines(method, -2, 2 + max(result.index), linestyles='--')
    plt.xlim(-2, 2 + max(result.index))
    plt.title('LOF局部離羣點檢測', fontsize=13)
    plt.ylabel('局部離羣因子', fontsize=15)
    plt.legend()
    plt.show()

def ensemble_lof(data, predict=None, k=5, groups=[], method='auto', weight=1):
    import pandas as pd
    # 判斷是否傳入測試數據,若沒有傳入則測試數據賦值爲訓練數據
    try:
        if predict == None:
            predict = data
    except Exception:
        pass
    data = pd.DataFrame(data);
    predict = pd.DataFrame(predict)
    # 數據標籤分組,默認獨立自成一組
    for i in range(data.shape[1]):
        if i not in pd.DataFrame(groups).values:
            groups += [[i]]
    # 擴充權值列表
    if type(weight) != list:
        weight = [weight]
        weight += [1] * (len(groups) - 1)
    else:
        weight += [1] * (len(groups) - len(weight))
    predict['local outlier factor'] = 0
    # 計算LOF離羣因子並根據特徵權重計算加權LOF得分
    for i in range(len(groups)):
        predict = localoutlierfactor(pd.DataFrame(data).iloc[:, groups[i]], predict, k, str(groups[i]))
        predict['local outlier factor'] += predict.iloc[:, -1] * weight[i]
    if method == 'auto':
        method = sum(weight)
    plot_lof(predict['local outlier factor'], method)
    # 根據離羣閾值劃分離羣點與正常點
    outliers = predict[predict['local outlier factor'] > method].sort_values(by='local outlier factor')
    inliers = predict[predict['local outlier factor'] <= method].sort_values(by='local outlier factor')
    return outliers, inliers
複製代碼

測試 5 仍然使用測試3的狀況進行分析,此時將經度、緯度設置爲獨立的特徵,分別對兩個維度數據進行識別(儘管單獨的緯度、經度數據彷佛沒有太大的實際意義)

import numpy as np
import pandas as pd

posi = pd.read_excel(r'已結束項目任務數據.xls')
lon = np.array(posi["任務gps經度"][:])  # 經度
lat = np.array(posi["任務gps 緯度"][:])  # 緯度
A = list(zip(lat, lon))  # 按照緯度-經度匹配

posi = pd.read_excel(r'會員信息數據.xlsx')
lon = np.array(posi["會員位置(GPS)經度"][:])  # 經度
lat = np.array(posi["會員位置(GPS)緯度"][:])  # 緯度
B = list(zip(lat, lon))  # 按照緯度-經度匹配

# 獲取會員對任務密度,取第5鄰域,閾值爲 100,權重分別爲5,1 
outliers5, inliers5 = ensemble_lof(A, B, k=5, method=100,weight = [5,1])

# 繪圖程序
plt.figure('LOF 異常分數加權模式')
plt.scatter(np.array(B)[:, 0], np.array(B)[:, 1], s=10, c='b', alpha=0.5)
plt.scatter(np.array(A)[:, 0], np.array(A)[:, 1], s=10, c='green', alpha=0.3)
plt.scatter(outliers5[0], outliers5[1], s=10 + outliers5['local outlier factor'], c='r', alpha=0.2)
plt.title('k = 5, method = 100')
複製代碼

比上面的各類識別模型更爲精確地識別出了可能的異常點

4.3.3 混合模式

混合模式適用於數據中有些特徵同等重要,有些特徵有重要性區別的狀況,即對 4.3.一、4.3.2 情形綜合進行考慮。同等重要的數據將使用投票表決模式,重要程度不一樣的數據使用加權模式並根據閾值轉換爲投票表決模式。程序上只需將兩部分混合使用便可,本文在此不作展現。

4.4 數據量太大,算法執行效率太低,有什麼改進方法嗎?

LOF 算法在檢測離羣點的過程當中,遍歷整個數據集以計算每一個點的 LOF 值,使得算法運算速度慢。同時,因爲數據正常點的數量通常遠遠多於離羣點的數量,而 LOF 方法經過比較全部數據點的 LOF 值判斷離羣程度,產生了大量不必的計算。所以,經過對原始數據進行修剪能夠有效提升 LOF 方法的計算效率。此外,實踐過程當中也發現經過數據集修剪後,能夠大幅度減小數據誤判爲離羣點的概率。這種基於聚類修剪得離羣點檢測方法稱爲 CLOF (Cluster-Based Local Outlier Factor) 算法。

基於 K-Means 的 CLOF 算法

在應用 LOF 算法前,先用 K-Means 聚類算法,將原始數據聚成 n\_clusters 簇。對其中的每一簇,計算簇的中心 C_i ,求出該簇中全部點到該中心的平均距離並記爲該簇的半徑 R_i 。對該類中全部點,若該點到簇中心的距離大於等於 R_i 則將其放入「離羣點候選集」 \bar{D} ,最後對 \bar{D} 中的數據使用 LOF 算法計算離羣因子。

設第 i 個簇中的點的個數爲 N(i) ,點集爲 \{X_{i1},X_{i2},\cdots,X_{i,N(i)}\}

中心和半徑的計算公式以下:

C_i=\frac{\displaystyle\sum^{N(i)}_{j=1}X_{ij}}{N(i)}\qquad R_i=\frac{\sqrt{\displaystyle\sum^{N(i)}_{j=1}(X_{ij}-C_i)^2}}{N(i)}\tag{20}

如何肯定最佳的 n\_clusters —— 肘部法則

K-Means算法經過指定聚類簇數n\_clusters及隨機生成聚類中心,對最靠近他們的對象進行迭代歸類,逐次更新各聚類中心的值,直到最好的聚類效果(代價函數值最小)。

對於n\_clusters的選取將直接影響算法的聚類效果。肘部法則將不一樣的n\_clusters值的成本函數值刻畫出來,隨着n\_clusters增大,每一個簇包含的樣本數會減小,樣本離其中心更接近,代價函數會減少。但隨着n\_clusters繼續增大,代價函數的改善程度不斷降低(在圖像中,代價函數曲線趨於平穩)。n\_clusters值增大過程當中,代價函數改善程度最大的位置對應的n\_clusters就是肘部,使用此n\_clusters通常能夠取得不錯的效果。但肘部法則的使用僅僅是從代價水平進行考慮,有時候還需結合實際考慮。

因爲離羣值樣本數量通常較少,若是聚類出來的簇中樣本量太少(如 1-4 個,但其餘簇有成百上千個樣本),則這種聚類簇不該進行修剪。

定義代價函數:

cost(n\_clusters)=\frac{\displaystyle\sum^{n\_clusters}_{i=1}\sqrt{\displaystyle\sum^{N(i)}_{j=1}(X_{ij}-C_i)^2}}{\mid\displaystyle\sum^{n\_clusters}_{i=1}N(i)\mid}
\tag{21}

CLOF局部離羣因子檢測算法流程

elbow_method(data,maxtest = 11)

  • 輸入:須要修剪的數據集,最大測試範圍(默認爲11)
  • 輸出:代價曲線圖
def elbow_method(data,maxtest = 11):
    '''使用肘部法則肯定$n_clusters$值'''
    from sklearn.cluster import KMeans
    from scipy.spatial.distance import cdist
    import numpy as np
    import matplotlib.pyplot as plt
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False
    ax = plt.figure(figsize=(8,4)).add_subplot(111)
    N_test = range(1, maxtest)
    # 代價函數值列表
    meandistortions = []
    for $n_clusters$ in N_test:
        model = KMeans($n_clusters$=$n_clusters$).fit(data)
        # 計算代價函數值
        meandistortions.append(sum(np.min(cdist(data, model.cluster_centers_, 'euclidean'), axis=1)) / len(data))
    plt.plot(N_test, meandistortions, 'bx-',alpha = 0.4)
    plt.xlabel('k')
    plt.ylabel('代價函數',fontsize= 12)
    plt.title('用肘部法則來肯定最佳的$n_clusters$值',fontsize= 12)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.set_xticks(np.arange(0, maxtest, 1))
    plt.show()
複製代碼

不一樣聚類簇數對應的代價曲線,圖中提示 2~4 可能爲理想值

kmeans(data, $n_clusters$, m):

  • 輸入:須要修剪的數據集,聚類簇數,最小修剪簇應包含的樣本點數
  • 輸出:修剪後的數據集
def kmeans(data, $n_clusters$, m):
    '''使用K-Means算法修剪數據集'''
    from sklearn.cluster import KMeans
    import numpy as np
    data_select = []
    model = KMeans($n_clusters$=$n_clusters$).fit(data)
    centeroids = model.cluster_centers_
    label_pred = model.labels_
    import collections
    for k, v in collections.Counter(label_pred).items():
        if v < m:
            data_select += np.array(data)[label_pred == k].tolist()
        else:
            distance = np.sqrt(((np.array(data)[label_pred == k] - centeroids[k]) ** 2).sum(axis=1))
            R = distance.mean()
            data_select += np.array(data)[label_pred == k][distance >= R].reshape(-1, np.array(data).shape[-1]).tolist()
    return np.array(data_select)
複製代碼

測試6 對B數據集進行修剪分析,B數據集共有 1877 條數據

# 肘部法則肯定最佳修剪集
elbow_method(B,maxtest = 11)
複製代碼

B數據集代價曲線:圖中提示 3 可能爲理想聚類簇數
# 根據上面設定的 $n_clusters$ 爲3,最小樣本量設置爲3
B_cut = kmeans(B, $n_clusters$ = 3, m = 3)
複製代碼

執行上述程序,原先包含 1877 條數據的 B 數據集修剪爲含有 719 條數據的較小的數據集。使用 LOF 算法進行離羣檢測,檢測結果以下:

# 獲取會員分佈密度,取第10鄰域,閾值爲3(LOF大於3認爲是離羣值)
outliers6, inliers6 = lof(B_cut, k=10, method=3)
# 繪圖程序
plt.figure('CLOF 離羣因子檢測')
plt.scatter(np.array(B)[:, 0], np.array(B)[:, 1], s=10, c='b', alpha=0.5)
plt.scatter(outliers6[0], outliers6[1], s=10 + outliers6['local outlier factor']*10, c='r', alpha=0.2)
plt.title('k = 10, method = 3')
複製代碼

使用 CLOF 檢測狀況,精度進一步獲得提高,減小了數據誤判率

修剪後的數據集 LOF 意義再也不那麼明顯,但離羣點的 LOF 仍然會是較大的值,而且 k 選取越大的值,判別效果越明顯。

4.5 如何提升識別精度

合理增大 k 值能顯著提升識別精度。但 k 值的增長會帶來沒必要要的計算、影響算法執行效率,也正所以本文所用的 k 值都取較小。合理選取 k 與閾值將是 LOF 成功與否的關鍵。

閾值選取 1.2,方便展現識別精度隨着k增大的提高過程,計算耗時也愈來愈長(僅取1次實驗,小樣本、小k值時效果不明顯,所以並無表出嚴格的遞增關係,但最後一圖能夠明顯看到耗時延長)。從氣泡大小也能夠看出,適當調高閾值也能提升識別精度,不必定依賴於k

5、後記

本文內容主要參考算法原文及筆者學習經驗進行總結。在異常識別領域,LOF 算法和 Isolation Forest 算法已經被指出是性能最優、識別效果最好的算法。對於經常使用的人羣密度(或其餘)的刻畫,LOF異常分數值不失爲一種「高端」方法(參考文獻[3]),相比傳統方法,更具備成熟的理論支撐。

後續有時間的話,筆者會根據此文方法,結合實際數據詳細進一步說明如何在數據處理中應用 LOF 算法。

6、參考文獻

[1] 維基百科. Moore–Penrose inverse[EB/OL]. [2018-6-2]. https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse

[2] Breunig M M, Kriegel H P, Ng R T, et al. LOF: identifying density-based local outliers[C]// ACM SIGMOD International Conference on Management of Data. ACM, 2000:93-104.

[3] 董天文,潘偉堤,戚銘珈,張曉敏. 「拍照賺錢」的任務訂價. 教育部中國大學生在線網站[J/OL]. [2017-11-13]. http://upload.univs.cn/2017/1113/1510570109509.pdf


做者:張柳彬 仇禮鴻

若有疑問,請聯繫QQ:965579168

轉載請聲明出處

相關文章
相關標籤/搜索