【免責聲明：用于教學資料整理】

目錄：

一. 用箱線圖檢測異常值

二. 使用局部異常因子法（LOF法）檢測異常值

三. 用聚類方法檢測異常值

四. 檢測時間序列數據中的異常值

五. 基于穩健馬氏距離檢測異常值

正文：

異常值，是指測量數據中的隨機錯誤或偏差，包括錯誤值或偏離均值的孤立點值。在數據處理中，異常值會極大的影響回歸或分類的效果。

為了避免異常值造成的損失，需要在數據預處理階段進行異常值檢測。另外，某些情況下，異常值檢測也可能是研究的目的，例如，數據造假的發現、電腦入侵的檢測等。

一、用箱線圖檢測異常值

在一條數軸上，以數據的上下四分位數（Q1-Q3）為界畫一個矩形盒子（中間50%的數據落在盒內）；在數據的中位數位置畫一條線段為中位線；用◇標記數據的均值；默認延長線不超過盒長的1.5倍，之外的點認為是異常值（用○標記）。

盒形圖的主要應用就是，剔除數據的異常值、判斷數據的偏態和尾重。

R語言實現，使用函數boxplot.stats()，基本格式為：

[stats, n, conf, out]=

boxplot.stats(x, coef=1.5, do.conf=TRUE, do.out=TRUE)

其中，x為數值向量（NA、NaN值將被忽略）；coef為盒須的長度為幾倍的IQR（盒長），默認為1.5；do.conf和do.out設置是否輸出conf和out

返回值：stats返回5個元素的向量值，包括盒須最小值、盒最小值、中位數、盒最大值、盒須最大值；n返回非缺失值的個數；conf返回中位數的95%置信區間；out返回異常值。

單變量異常值檢測：

set.seed(2016)

x<-rnorm(100) #生成100個服從N(0,1)的隨機數

summary(x) #x的匯總信息

Min. 1st Qu. Median Mean 3rd Qu. Max.

-2.7910 -0.7173 -0.2662 -0.1131 0.5917 2.1940

boxplot.stats(x) #用箱線圖檢測x中的異常值

$stats

[1] -2.5153136 -0.7326879 -0.2662071 0.5929206 2.1942200

$n

[1] 100

$conf

[1] -0.47565320 -0.05676092

$out

[1] -2.791471

boxplot(x) #繪制箱線圖

多變量異常值檢測：

x<-rnorm(100)

y<-rnorm(100)

df<-data.frame(x,y) #用x,y生成兩列的數據框

head(df)

x y

1 0.41452353 0.4852268

2 -0.47471847 0.6967688

3 0.06599349 0.1855139

4 -0.50247778 0.7007335

5 -0.82599859 0.3116810

6 0.16698928 0.7604624

#尋找x為異常值的坐標位置

a<-which(x %in% boxplot.stats(x)$out)

a

[1] 78 81 92

#尋找y為異常值的坐標位置

b<-which(y %in% boxplot.stats(y)$out)

b

[1] 27 37

intersect(a,b) #尋找變量x,y都為異常值的坐標位置

integer(0)

plot(df) #繪制x, y的散點圖

p2<-union(a,b) #尋找變量x或y為異常值的坐標位置

[1] 78 81 92 27 37

points(df[p2,],col="red",pch="x",cex=2) #標記異常值

二、使用局部異常因子法（LOF法）檢測異常值

局部異常因子法（LOF法），是一種基于概率密度函數識別異常值的算法。LOF算法只對數值型數據有效。

算法原理：將一個點的局部密度與其周圍的點的密度相比較，若前者明顯的比后者小（LOF值大于1），則該點相對于周圍的點來說就處于一個相對比較稀疏的區域，這就表明該點是一個異常值。

R語言實現：使用DMwR或dprep包中的函數lofactor()，基本格式為：

lofactor(data, k)

其中，data為數值型數據集；k為用于計算局部異常因子的鄰居數量。

library(DMwR)

iris2<-iris[,1:4] #只選數值型的前4列

head(iris2)

Sepal.Length Sepal.Width Petal.Length Petal.Width

1 5.1 3.5 1.4 0.2

2 4.9 3.0 1.4 0.2

3 4.7 3.2 1.3 0.2

4 4.6 3.1 1.5 0.2

5 5.0 3.6 1.4 0.2

6 5.4 3.9 1.7 0.4

out.scores<-lofactor(iris2,k=10) #計算每個樣本的LOF值

plot(density(out.scores)) #繪制LOF值的概率密度圖

#LOF值排前5的數據作為異常值，提取其樣本號

out<-order(out.scores,decreasing=TRUE)[1:5]

out

[1] 42 107 23 16 99

iris2[out,] #異常值數據

Sepal.Length Sepal.Width Petal.Length Petal.Width

42 4.5 2.3 1.3 0.3

107 4.9 2.5 4.5 1.7

23 4.6 3.6 1.0 0.2

16 5.7 4.4 1.5 0.4

99 5.1 2.5 3.0 1.1

對鳶尾花數據進行主成分分析，并利用產生的前兩個主成分繪制成雙標圖來顯示異常值：

n<-nrow(iris2) #樣本數

n

[1] 150

labels<-1:n #用數字1-n標注

labels[-out]<-"." #非異常值用"."標注

biplot(prcomp(iris2),cex=0.8,xlabs=labels)

說明：函數prcomp()對數據集iris2做主成份分析，biplot()取主成份分析結果的前兩列數據即前兩個主成份繪制雙標圖。上圖中，x軸和y軸分別代表第一、二主成份，箭頭指向了原始變量名，其中5個異常值分別用對應的行號標注。

也可以通過函數pairs()繪制散點圖矩陣來顯示異常值，其中異常值用紅色的"+"標注：

pchs<-rep(".",n)

pchs[out]="+"

cols<-rep("black",n)

cols[out]<-"red"

pairs(iris2,pch=pchs,col=cols)

注：另外，Rlof包中函數lof()可實現相同的功能，并且支持并行計算和選擇不同距離。

三、用聚類方法檢測異常值

通過把數據聚成類，將那些不屬于任何一類的數據作為異常值。比如，使用基于密度的聚類DBSCAN，如果對象在稠密區域緊密相連，則被分組到一類；那些不會被分到任何一類的對象就是異常值。

也可以用k-means算法來檢測異常值：將數據分成k組，通過把它們分配到最近的聚類中心。然后，計算每個對象到聚類中心的距離（或相似性），并選擇最大的距離作為異常值。

kmeans.result<-kmeans(iris2,centers=3) #kmeans聚類為3類

kmeans.result$centers #輸出聚類中心

Sepal.Length Sepal.Width Petal.Length Petal.Width

1 5.901613 2.748387 4.393548 1.433871

2 5.006000 3.428000 1.462000 0.246000

3 6.850000 3.073684 5.742105 2.071053

kmeans.result$cluster #輸出聚類結果

[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

[30] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 3 1 1 1 1 1

[59] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1

[88] 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 3 3 3 3 1 3 3 3 3 3 3 1 1 3

[117] 3 3 3 1 3 1 3 1 3 3 1 1 3 3 3 3 3 1 3 3 3 3 1 3 3 3 1 3 3

[146] 3 1 3 3 1

#centers返回每個樣本對應的聚類中心樣本

centers <- kmeans.result$centers[kmeans.result$cluster, ]

#計算每個樣本到其聚類中心的距離

distances<-sqrt(rowSums((iris2-centers)^2))

#找到距離最大的5個樣本，認為是異常值

out<-order(distances,decreasing=TRUE)[1:5]

out #異常值的樣本號

[1] 99 58 94 61 119

iris2[out,] #異常值

Sepal.Length Sepal.Width Petal.Length Petal.Width

99 5.1 2.5 3.0 1.1

58 4.9 2.4 3.3 1.0

94 5.0 2.3 3.3 1.0

61 5.0 2.0 3.5 1.0

119 7.7 2.6 6.9 2.3

#繪制聚類結果

plot(iris2[,c("Sepal.Length","Sepal.Width")],pch="o",col=kmeans.result$cluster,cex=0.3)

#聚類中心用"*"標記

points(kmeans.result$centers[,c("Sepal.Length", "Sepal.Width")], col=1:3, pch=8, cex=1.5)

#異常值用"+"標記

points(iris2[out,c("Sepal.Length", "Sepal.Width")], pch="+", col=4, cex=1.5)

四、檢測時間序列數據中的異常值

對時間序列數據進行異常值檢測，先用函數stl()進行穩健回歸分解，再識別異常值。

函數stl()，基于局部加權回歸散點平滑法(LOESS)，對時間序列數據做穩健回歸分解，分解為季節性、趨勢性、不規則性三部分。

f<-stl(AirPassengers,"periodic",robust=TRUE)

#weights返回穩健性權重，以控制數據中異常值產生的影響

out<-which(f$weights < 1e-8) #找到異常值

out

[1] 79 91 92 102 103 104 114 115 116 126 127 128 138 139 140

#設置繪圖布局的參數

op<-par(mar=c(0,4,0,3), oma=c(5,0,4,0), mfcol=c(4,1))

plot(f,set.pars=NULL)

#time.series返回分解為三部分的時間序列

> head(f$time.series,3)

seasonal trend remainder

[1,] -16.519819 123.1857 5.3341624

[2,] -27.337882 123.4214 21.9164399

[3,] 9.009778 123.6572 -0.6670047

sts<-f$time.series

#用紅色"x"標記異常值

points(time(sts)[out], 0.8*sts[,"remainder"][out], pch="x", col="red")

par(op)

五、基于穩健馬氏距離檢測異常值

檢驗異常值的基本思路是觀察各樣本點到樣本中心的距離，若某些樣本點的距離太大，就可以判斷是異常值。

若使用歐氏距離，則具有明顯的缺點：將樣本不同屬性（即各指標變量）之間的差別等同看待。而馬氏距離則不受量綱的影響，并且在多元條件下，還考慮到了變量之間的相關性。

對均值為μ，協方差矩陣為Σ的多變量向量，其馬氏距離為

(x-μ)TΣ-1(x-μ)

但是傳統的馬氏距離檢測方法是不穩定的，因為個別異常值會把均值向量和協方差矩陣向自己方向吸引，這就導致馬氏距離起不了檢測異常值的所用。解決方法是利用迭代思想構造一個穩健的均值和協方差矩陣估計量，然后計算穩健馬氏距離，這樣異常值就能正確地被識別出來。

用mvoutlier包實現，

library(mvoutlier)

set.seed(2016)

x<-cbind(rnorm(80),rnorm(80))

y<-cbind(rnorm(10,5,1), rnorm(10,5,1)) #噪聲數據

z<-rbind(x,y)

res1<-uni.plot(z) #一維數據的異常值檢驗

#返回outliers標記各樣本是否為異常值，md返回數據的穩健馬氏距離

which(res1$outliers==TRUE) #返回異常值的樣本號

[1] 81 82 83 84 85 86 87 88 89 90

res2<-aq.plot(z) #基于穩健馬氏距離的多元異常值檢驗

which(res2$outliers==TRUE) #返回異常值的樣本號

[1] 81 82 83 84 85 86 87 88 89 90

上圖為在一維空間中觀察樣本數據。

說明：圖1-1為原始數據；圖1-2的X軸為各樣本的穩健馬氏距離排序，Y軸為距離的經驗分布，紅色曲線為卡方分布，藍色垂線表示閥值，在閥值右側的樣本判斷為異常值；圖2-1和2-2均是用不同顏色來表示異常值，只是閥值略有不同。

若數據的維數過高，則上述距離不再有很大意義（例如基因數據有幾千個變量，數據之間變得稀疏）。此時可以融合主成份降維的思路來進行異常值檢驗。mvoutlier包中提供了函數pcout()來對高維數據進行異常值檢驗。

data(swiss) #使用swiss數據集

res3<-pcout(swiss)

#返回wfinal01標記是否為異常值，0表示是

which(res3$wfinal01==0) #返回異常值的樣本號

Delemont Franches-Mnt Porrentruy Broye

2 3 6 7

Glane Gruyere Sarine Veveyse

8 9 10 11

La Vallee Conthey Entremont Herens

19 31 32 33

Martigwy Monthey St Maurice Sierre

34 35 36 37

Sion V. De Geneve

38 45

注：對于分類數據，一個快速穩定的異常檢測的策略是AVF (Attribute Value Frequency)算法。

主要參考文獻：

《R語言-異常值處理1-3》，銀河統計學，博客園

http://www.cnblogs.com/cloudtj/category/780800.html

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的剔除异常值栅格计算器_R语言系列 数据清洗3 异常值处理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。