R语言学习系列16-异常值处理

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.LengthSepal.WidthPetal.LengthPetal.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(\ol=kmeans.result$cluster,cex=0.3)

#聚类中心用\标记

points(kmeans.result$centers[,c(\Width\

#异常值用\标记

points(iris2[out,c(\=\

四、检测时间序列数据中的异常值

对时间序列数据进行异常值检测，先用函数stl()进行稳健回归分解，再识别异常值。

函数stl()，基于局部加权回归散点平滑法(LOESS)，对时间序列数据做稳健回归分解，分解为季节性、趋势性、不规则性三部分。

f<-stl(AirPassengers,\#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

#用红色\标记异常值

points(time(sts)[out], 0.8*sts[,\\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

上图为在一维空间中观察样本数据。