【免責(zé)聲明：用于教學(xué)資料整理】

目錄：

一. 用箱線圖檢測異常值

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

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

四. 檢測時(shí)間序列數(shù)據(jù)中的異常值

五. 基于穩(wěn)健馬氏距離檢測異常值

正文：

異常值，是指測量數(shù)據(jù)中的隨機(jī)錯(cuò)誤或偏差，包括錯(cuò)誤值或偏離均值的孤立點(diǎn)值。在數(shù)據(jù)處理中，異常值會極大的影響回歸或分類的效果。

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

一、用箱線圖檢測異常值

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

盒形圖的主要應(yīng)用就是，剔除數(shù)據(jù)的異常值、判斷數(shù)據(jù)的偏態(tài)和尾重。

R語言實(shí)現(xiàn)，使用函數(shù)boxplot.stats()，基本格式為：

[stats, n, conf, out]=

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

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

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

單變量異常值檢測：

set.seed(2016)

x<-rnorm(100) #生成100個(gè)服從N(0,1)的隨機(jī)數(shù)

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生成兩列的數(shù)據(jù)框

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為異常值的坐標(biāo)位置

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

a

[1] 78 81 92

#尋找y為異常值的坐標(biāo)位置

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

b

[1] 27 37

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

integer(0)

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

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

[1] 78 81 92 27 37

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

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

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

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

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

lofactor(data, k)

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

library(DMwR)

iris2<-iris[,1:4] #只選數(shù)值型的前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) #計(jì)算每個(gè)樣本的LOF值

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

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

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

out

[1] 42 107 23 16 99

iris2[out,] #異常值數(shù)據(jù)

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

對鳶尾花數(shù)據(jù)進(jìn)行主成分分析，并利用產(chǎn)生的前兩個(gè)主成分繪制成雙標(biāo)圖來顯示異常值：

n<-nrow(iris2) #樣本數(shù)

n

[1] 150

labels<-1:n #用數(shù)字1-n標(biāo)注

labels[-out]<-"." #非異常值用"."標(biāo)注

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

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

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

pchs<-rep(".",n)

pchs[out]="+"

cols<-rep("black",n)

cols[out]<-"red"

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

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

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

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

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

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 #輸出聚類結(jié)果

[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返回每個(gè)樣本對應(yīng)的聚類中心樣本

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

#計(jì)算每個(gè)樣本到其聚類中心的距離

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

#找到距離最大的5個(gè)樣本，認(rèn)為是異常值

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

#繪制聚類結(jié)果

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

#聚類中心用"*"標(biāo)記

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

#異常值用"+"標(biāo)記

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

四、檢測時(shí)間序列數(shù)據(jù)中的異常值

對時(shí)間序列數(shù)據(jù)進(jìn)行異常值檢測，先用函數(shù)stl()進(jìn)行穩(wěn)健回歸分解，再識別異常值。

函數(shù)stl()，基于局部加權(quán)回歸散點(diǎn)平滑法(LOESS)，對時(shí)間序列數(shù)據(jù)做穩(wěn)健回歸分解，分解為季節(jié)性、趨勢性、不規(guī)則性三部分。

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

#weights返回穩(wěn)健性權(quán)重，以控制數(shù)據(jù)中異常值產(chǎn)生的影響

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

out

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

#設(shè)置繪圖布局的參數(shù)

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返回分解為三部分的時(shí)間序列

> 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"標(biāo)記異常值

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

par(op)

五、基于穩(wěn)健馬氏距離檢測異常值

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

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

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

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

但是傳統(tǒng)的馬氏距離檢測方法是不穩(wěn)定的，因?yàn)閭€(gè)別異常值會把均值向量和協(xié)方差矩陣向自己方向吸引，這就導(dǎo)致馬氏距離起不了檢測異常值的所用。解決方法是利用迭代思想構(gòu)造一個(gè)穩(wěn)健的均值和協(xié)方差矩陣估計(jì)量，然后計(jì)算穩(wěn)健馬氏距離，這樣異常值就能正確地被識別出來。

用mvoutlier包實(shí)現(xiàn)，

library(mvoutlier)

set.seed(2016)

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

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

z<-rbind(x,y)

res1<-uni.plot(z) #一維數(shù)據(jù)的異常值檢驗(yàn)

#返回outliers標(biāo)記各樣本是否為異常值，md返回?cái)?shù)據(jù)的穩(wěn)健馬氏距離

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

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

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

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

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

上圖為在一維空間中觀察樣本數(shù)據(jù)。

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

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

data(swiss) #使用swiss數(shù)據(jù)集

res3<-pcout(swiss)

#返回wfinal01標(biāo)記是否為異常值，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

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

主要參考文獻(xiàn)：

《R語言-異常值處理1-3》，銀河統(tǒng)計(jì)學(xué)，博客園

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

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。