搜索
16.异常值处理目录：一.用箱线图检测异常值二.使用局部异常因子法（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.1stQu.MedianMean3rdQu.Max.-2.7910-0.7173-0.2662-0.11310.59172.1940boxplot.stats(x)#用箱线图检测x中的异常值$stats[1]-2.5153136-0.7326879-0.26620710.59292062.1942200$n[1]100$conf[1]-0.47565320-0.05676092$out[1]-2.791471boxplot(x)#绘制箱线图多变量异常值检测：x-rnorm(100)y-rnorm(100)df-data.frame(x,y)#用x,y生成两列的数据框head(df)xy10.414523530.48522682-0.474718470.696768830.065993490.18551394-0.502477780.70073355-0.825998590.311681060.166989280.7604624#寻找x为异常值的坐标位置a-which(x%in%boxplot.stats(x)$out)a[1]788192#寻找y为异常值的坐标位置b-which(y%in%boxplot.stats(y)$out)b[1]2737intersect(a,b)#寻找变量x,y都为异常值的坐标位置integer(0)plot(df)#绘制x,y的散点图p2-union(a,b)#寻找变量x或y为异常值的坐标位置[1]7881922737points(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.LengthSepal.WidthPetal.LengthPetal.Width15.13.51.40.224.93.01.40.234.73.21.30.244.63.11.50.255.03.61.40.265.43.91.70.4out.scores-lofactor(iris2,k=10)#计算每个样本的LOF值plot(density(out.scores))#绘制LOF值的概率密度图#LOF值排前5的数据作为异常值，提取其样本号out-order(out.scores,decreasing=TRUE)[1:5]out[1]42107231699iris2[out,]#异常值数据Sepal.LengthSepal.WidthPetal.LengthPetal.Width424.52.31.30.31074.92.54.51.7234.63.61.00.2165.74.41.50.4995.12.53.01.1对鸢尾花数据进行主成分分析，并利用产生的前两个主成分绘制成双标图来显示异常值：n-nrow(iris2)#样本数n[1]150labels-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]-redpairs(iris2,pch=pchs,col=cols)注：另外，Rlof包中函数lof()可实现相同的功能，并且支持并行计算和选择不同距离。三、用聚类方法检测异常值通过把数据聚成类，将那些不属于任何一类的数据作为异常值。比如，使用基于密度的聚类DBSCAN，如果对象在稠密区域紧密相连，则被分组到一类；那些不会被分到任何一类的对象就是异常值。也可以用k-means算法来检测异常值：将数据分成k组，通过把它们分配到最近的聚类中心。然后，计算每个对象到聚类中心的距离（或相似性），并选择最大的距离作为异常值。kmeans.result-kmeans(iris2,centers=3)#kmeans聚类为3类kmeans.result$centers#输出聚类中心Sepal.LengthSepal.WidthPetal.LengthPetal.Width15.9016132.7483874.3935481.43387125.0060003.4280001.4620000.24600036.8500003.0736845.7421052.071053kmeans.result$cluster#输出聚类结果[1]22222222222222222222222222222[30]22222222222222222222211311111[59]11111111111111111113111111111[88]11111111111113133331333333113[117]33313131331133333133331333133[146]31331#centers返回每个样本对应的聚类中心样本centers-kmeans.result$centers[kmeans.result$cluster,]#计算每个样本到其聚类中心的距离distances-sqrt(rowSums((iris2-centers)^2))#找到距离最大的5个样本，认为是异常值out-order(distances,decreasing=TRUE)[1:5]out#异常值的样本号[1]99589461119iris2[out,]#异常值Sepal.LengthSepal.WidthPetal.LengthPetal.Width995.12.53.01.1584.92.43.31.0945.02.33.31.0615.02.03.51.01197.72.66.92.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$weights1e-8)#找到异常值out[1]799192102103104114115116126127128138139140#设置绘图布局的参数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)seasonaltrendremainder[1,]-16.519819123.18575.3341624[2,]-27.337882123.421421.9164399[3,]9.009778123.6572-0.6670047sts-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]81828384858687888990res2-aq.plot(z)#基于稳健马氏距离的多元异常值检验which(res2$outliers==TRUE)#