R语言学习系列16异常值处理Word下载.docx

1、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.1940boxplot.stats（x） #用箱线图检测x中的异常值$stats1 -2.5153136 -0.7326879 -0.2662071 0.5929206 2.1942200$n1 100$conf1 -0.

2、47565320 -0.05676092$out1 -2.791471boxplot（x） #绘制箱线图多变量异常值检测：-rnorm（100）ydf-data.frame（x,y） #用x,y生成两列的数据框head（df） x y1 0.41452353 0.48522682 -0.47471847 0.69676883 0.06599349 0.18551394 -0.50247778 0.70073355 -0.82599859 0.31168106 0.16698928 0.7604624#寻找x为异常值的坐标位置a-which（x %in% boxplot.stats（x）$out

3、） a1 78 81 92#寻找y为异常值的坐标位置b-which（y %in% boxplot.stats（y）$out）b1 27 37intersect（a,b） #寻找变量x,y都为异常值的坐标位置integer（0）plot（df） #绘制x, y的散点图p2-union（a,b） #寻找变量x或y为异常值的坐标位置1 78 81 92 27 37points（dfp2,col=red,pch=x,cex=2） #标记异常值二、使用局部异常因子法（LOF法）检测异常值局部异常因子法（LOF法），是一种基于概率密度函数识别异常值的算法。LOF算法只对数值型数据有效。算法原理：将一个点的

4、局部密度与其周围的点的密度相比较，若前者明显的比后者小（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.Width1 5.1 3.5 1.4 0.22 4.9 3.0 1.4 0.2

5、3 4.7 3.2 1.3 0.24 4.6 3.1 1.5 0.25 5.0 3.6 1.4 0.26 5.4 3.9 1.7 0.4out.scores-lofactor（iris2,k=10） #计算每个样本的LOF值plot（density（out.scores） #绘制LOF值的概率密度图#LOF值排前5的数据作为异常值，提取其样本号out-order（out.scores,decreasing=TRUE）1:5out1 42 107 23 16 99iris2out, #异常值数据42 4.5 2.3 1.3 0.3107 4.9 2.5 4.5 1.723 4.6 3.6 1.0

6、 0.216 5.7 4.4 1.5 0.499 5.1 2.5 3.0 1.1对鸢尾花数据进行主成分分析，并利用产生的前两个主成分绘制成双标图来显示异常值：n-nrow（iris2） #样本数n1 150labels-1:n #用数字1-n标注labels-out-. #非异常值用标注biplot（prcomp（iris2）,cex=0.8,xlabs=labels） 说明：函数prcomp（）对数据集iris2做主成份分析，biplot（）取主成份分析结果的前两列数据即前两个主成份绘制双标图。上图中，x轴和y轴分别代表第一、二主成份，箭头指向了原始变量名，其中5个异常值分别用对应的行号标注

7、。 也可以通过函数pairs（）绘制散点图矩阵来显示异常值，其中异常值用红色的+标注：pchs-rep（,n）pchsout=colsblackcolsoutpairs（iris2,pch=pchs,col=cols） 注：另外，Rlof包中函数lof（）可实现相同的功能，并且支持并行计算和选择不同距离。三、用聚类方法检测异常值 通过把数据聚成类，将那些不属于任何一类的数据作为异常值。比如，使用基于密度的聚类DBSCAN，如果对象在稠密区域紧密相连，则被分组到一类；那些不会被分到任何一类的对象就是异常值。也可以用k-means算法来检测异常值：将数据分成k组，通过把它们分配到最近的聚类中心。然

8、后，计算每个对象到聚类中心的距离（或相似性），并选择最大的距离作为异常值。kmeans.result-kmeans（iris2,centers=3） #kmeans聚类为3类kmeans.result$centers #输出聚类中心1 5.901613 2.748387 4.393548 1.4338712 5.006000 3.428000 1.462000 0.2460003 6.850000 3.073684 5.742105 2.071053kmeans.result$cluster #输出聚类结果 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

9、 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 3117 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 3146 3 1 3 3 1#centers返回每个样本对应的聚类中心样本cent

10、ers - kmeans.result$centerskmeans.result$cluster, #计算每个样本到其聚类中心的距离distances-sqrt（rowSums（iris2-centers）2）#找到距离最大的5个样本，认为是异常值-order（distances,decreasing=TRUE）1:out #异常值的样本号1 99 58 94 61 119iris2out, #异常值58 4.9 2.4 3.3 1.094 5.0 2.3 3.3 1.061 5.0 2.0 3.5 1.0119 7.7 2.6 6.9 2.3#绘制聚类结果plot（iris2,c（Sepal

11、.Length,Sepal.Width）,pch=o,col=kmeans.result$cluster,cex=0.3）#聚类中心用*标记points（kmeans.result$centers,c（, ）, col=1:3, pch=8, cex=1.5）#异常值用points（iris2out,c（）, pch=, col=4, cex=1.5）四、检测时间序列数据中的异常值对时间序列数据进行异常值检测，先用函数stl（）进行稳健回归分解，再识别异常值。函数stl（），基于局部加权回归散点平滑法（LOESS），对时间序列数据做稳健回归分解，分解为季节性、趋势性、不规则性三部分。f-stl

12、（AirPassengers,periodic,robust=TRUE）#weights返回稳健性权重，以控制数据中异常值产生的影响-which（f$weights 1e-8） #找到异常值1 79 91 92 102 103 104 114 115 116 126 127 128 138 139 140#设置绘图布局的参数op head（f$time.series,3） seasonal trend remainder1, -16.519819 123.1857 5.33416242, -27.337882 123.4214 21.91643993, 9.009778 123.6572 -0

13、.6670047sts-f$time.series#用红色标记异常值points（time（sts）out, 0.8*sts,remainderout, pch=, col=）par（op）五、基于稳健马氏距离检测异常值 检验异常值的基本思路是观察各样本点到样本中心的距离，若某些样本点的距离太大，就可以判断是异常值。若使用欧氏距离，则具有明显的缺点：将样本不同属性（即各指标变量）之间的差别等同看待。而马氏距离则不受量纲的影响，并且在多元条件下，还考虑到了变量之间的相关性。对均值为，协方差矩阵为的多变量向量，其马氏距离为（x-）T-1（x-）但是传统的马氏距离检测方法是不稳定的，因为个别异常值会

14、把均值向量和协方差矩阵向自己方向吸引，这就导致马氏距离起不了检测异常值的所用。解决方法是利用迭代思想构造一个稳健的均值和协方差矩阵估计量，然后计算稳健马氏距离，这样异常值就能正确地被识别出来。 用mvoutlier包实现，library（mvoutlier）-cbind（rnorm（80）,rnorm（80）-cbind（rnorm（10,5,1）, rnorm（10,5,1） #噪声数据z-rbind（x,y）res1-uni.plot（z） #一维数据的异常值检验#返回outliers标记各样本是否为异常值，md返回数据的稳健马氏距离which（res1$outliers=TRUE） #返

15、回异常值的样本号 1 81 82 83 84 85 86 87 88 89 90res2-aq.plot（z） #基于稳健马氏距离的多元异常值检验which（res2$outliers=TRUE） #返回异常值的样本号上图为在一维空间中观察样本数据。说明：图1-1为原始数据；图1-2的X轴为各样本的稳健马氏距离排序，Y轴为距离的经验分布，红色曲线为卡方分布，蓝色垂线表示阀值，在阀值右侧的样本判断为异常值；图2-1和2-2均是用不同颜色来表示异常值，只是阀值略有不同。若数据的维数过高，则上述距离不再有很大意义（例如基因数据有几千个变量，数据之间变得稀疏）。此时可以融合主成份降维的思路来进行异常值

16、检验。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，银河统计学，博客园