R语言--异常值检测

自编函数,boxplot()原理

outlier.IQR <- function(x, multiple = 1.5, replace = FALSE, revalue = NA) { 
  q <- quantile(x, na.rm = TRUE) #四分位间距3倍间距以外的认为是离群值
  IQR <- q[4] - q[2]
  x1 <- which(x < q[2] - multiple * IQR | x > q[4] + multiple * IQR)
  x2 <- x[x1]
  if (length(x2) > 0) outlier <- data.frame(location = x1, value = x2)
  else outlier <- data.frame(location = 0, value = 0)
  if (replace == TRUE) {
    x[x1] <- revalue
  }
  return(list(new.value = x, outlier = outlier))
}

结果输出为列表，分别为 outlier.IQR()$new.value 和 outlier.IQR()$outlier。前者为异常值替换后的新向量，后者为原向量中异常值及其所在位置。

异常检测，主要内容如下：

（1）单变量的异常检测

（2）使用LOF（local outlier factor，局部异常因子）进行异常检测

（3）通过聚类进行异常检测

（4）对时间序列进行异常检测

单变量异常检测

本部分展示了一个单变量异常检测的例子，并且演示了如何将这种方法应用在多元数据上。在该例中，单变量异常检测通过boxplot.stats()函数实现，并且返回产生箱线图的统计量。在返回的结果中，有一个部分是out，它结出了异常值的列表。更明确点，它列出了位于极值之外的胡须。参数coef可以控制胡须延伸到箱线图外的远近。在R中，运行?boxplot.stats可获取更详细的信息。

如图呈现了一个箱线图，其中有四个圈是异常值。

> set.seed(1234)
> x <- rnorm(1000)
> summary(x)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-3.39606 -0.67325 -0.03979 -0.02660  0.61582  3.19590 
> boxplot.stats(x)$out
 [1]  3.043766 -2.732220 -2.855759  2.919140 -3.233152 -2.651741
 [7] -3.396064  3.195901 -2.729680 -2.704203 -2.864347 -2.661346
[13]  2.705775 -2.906674 -2.874042 -2.757050 -2.739754
> y=rep(1,1000)
> z = data.frame(x,y)
> g <- ggplot(z,aes(y=x,x=y))
> g+geom_boxplot()

如上的单变量异常检测可以用来发现多元数据中的异常值，通过简单搭配的方式。在下例中，我们首先产生一个数据框df，它有两列x和y。之后，异常值分别从x和y检测出来。然后，我们获取两列都是异常值的数据作为异常数据。

在下图中，异常值用红色标记为”+”

> y = rnorm(1000)
> df <- data.frame(x,y)
> rm(x,y)
> head(df)
           x          y
1 -1.2070657 -1.2053334
2  0.2774292  0.3014667
3  1.0844412 -1.5391452
4 -2.3456977  0.6353707
5  0.4291247  0.7029518
6  0.5060559 -1.9058829

> attach(df)
> #find the index of outliers from x
> (a<- which(x  %in% boxplot.stats(x)$out))
 [1] 178 181 192 227 237 382 392 486 487 517 558 717 771 788 901 949
[17] 967
> 
> #find the index of outliers from y
> 
> (b <- which(y %in% boxplot.stats(y)$out))
[1] 121 233 317 359 517 660 815
> 

> detach(df)
> #outliers in both x and y
> (outlier.list1 <- intersect(a,b))
 [1] 517 
> 

plot(df)
points(df[outlier.list1,],col="red",cex=2.5,pch="+")

#或者用ggplot2
z = vector()
for(i in 1:1000){
  z[i]= ifelse (is.element(i,outlier.list1) ,1,2)
}
tt = cbind(df,z)
f<- ggplot(tt,aes(x = x,y=y))
f+geom_point(color=z,alpha=0.5)

类似的，我们也可以将x或y为异常值的数据标记为异常值。下图，异常值用’x’标记为蓝色。

#outliers in either x or y
(outlier.list2<- union(a,b))
plot(df)
points(df[outlier.list2,],col="blue",pch="x",cex=2)

(outlier.list1 <- union(a,b))
z = vector()
for(i in 1:1000){
  z[i]= ifelse (is.element(i,outlier.list1) ,1,2)
}
tt = cbind(df,z)

f<- ggplot(tt,aes(x = x,y=y))
f+geom_point(color=z,alpha=0.5)

当有三个以上的变量时，最终的异常值需要考虑单变量异常检测结果的多数表决。当选择最佳方式在真实应用中进行搭配时，需要涉及领域知识。

使用LOF（local outlier factor，局部异常因子）进行异常检测

LOF（局部异常因子）是用于识别基于密度的局部异常值的算法。使用LOF，一个点的部密度会与它的邻居进行比较。如果前者明显低于后者（有一个大于1 的LOF值），该位于一个稀疏区域，对于它的邻居而言，这就表明，该点是一个异常值。LOF的缺点就是它只对数值数据有效。

lofactor()函数使用LOF算法计算局部异常因子，
并且它在DMwR和dprep包中是可用的。下面将介绍一个使用LOF进行异常检测的例子，k是用于计算局部异常因子的邻居数量。下图呈现了一个异常值得分的密度图。

library(DMwR)
#remove"Species",which is a categorical column
iris2 <- iris[,1:4]
outlier.scores <- lofactor(iris2,k=5)
plot(density(outlier.scores))

#或者ggplot2
ggplot(as.data.frame(outlier.scores),aes(x=outlier.scores))+geom_density()

#pick top 5 as outliers
> outliers <- order(outlier.scores,decreasing = T)[1:5]
> #who are outliers
> print(outliers)
[1]  42 107  23 110  63
> 
> print(iris2[outliers,])
    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
110          7.2         3.6          6.1         2.5
63           6.0         2.2          4.0         1.0
> 

接着，我们结合前两个主成分的双标图呈现异常值

n <- nrow(iris2)
labels <- 1:n
labels[-outliers] <- "."
biplot(prcomp(iris2),cex = 0.8,xlabs = labels)

在如上代码中，prcomp()执行了一个主成分分析，并且biplot()使用前两个主成分画出了这些数据。在上图中，x和y轴分别代表第一和第二个主成份，箭头表示了变量，5个异常值用它们的行号标记出来了。

我们也可以如下使用pairsPlot显示异常值，这里的异常值
用”+”标记为红色。

pch<- rep(".",n)
pch[outliers]<- "+"
col <- rep("black",n)
col[outliers] <- "red"
pairs(iris2,pch = pch,col = col)

Rlof包，对LOF算法的并行实现。它的用法与lofactor()相似，但是lof()有两个附加的性，即支持k的多元值和距离度量的几种选择。如下是lof()的一个例子。在计算异常值得分后，异常值可以通过选择前几个检测出来。注意，目前包Rlof的版本在MacOS X和Linux环境下工作，但并不在windows环境下工作，因为它要依赖multicore包用于并行计算。

library(Rlof)
outlier.scores <- lof(iris2,k=5)
#try with different number of neighbors(k=5,6,7,8,9 and 10)
outlier.scores <- lof(iris2,k=c(5:10))

通过聚类进行异常检测

另外一种异常检测的方法是聚类。通过把数据聚成类，将那些不属于任务一类的数据
作为异常值。比如，使用基于密度的聚类DBSCAN，如果对象在稠密区域紧密相连，它们将被分组到一类。因此，那些不会被分到任何一类的对象就是异常值。

我们也可以使用k-means算法来检测异常。使用k-means算法，数据被分成k组，通过把它们分配到最近的聚类中心。然后，我们能够计算每个对象到聚类中心的距离（或相似性），并且选择最大的距离作为异常值。

如下是一个基于k-means算法在iris数据上实现在异常检测。

#remove species from the data to cluster
> iris2 <- iris[,1:4]
> kmeans.result <- kmeans(iris2,centers=3)
> #cluster centers
> kmeans.result$centers
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1     5.006000    3.428000     1.462000    0.246000
2     5.901613    2.748387     4.393548    1.433871
3     6.850000    3.073684     5.742105    2.071053
> 
> #cluster IDs
> kmeans.result$cluster
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 3 2 2 2
 [57] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2
 [85] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3 3 2 3 3 3 3 3
[113] 3 2 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 3 3 3 3 2 3
[141] 3 3 2 3 3 3 2 3 3 2
> 

#calculate distance between objects and cluster centers
> centers <- kmeans.result$centers[kmeans.result$cluster,]
> distances <- kmeans.result$centers[kmeans.result$cluster,]
> #pick top 5 largest distances
> outliers <- order(distances,decreasing = T)[1:5]
> #who are outliers
> print(outliers)
[1]  53  78 101 103 104
> print(iris2[outliers,])
    Sepal.Length Sepal.Width Petal.Length Petal.Width
53           6.9         3.1          4.9         1.5
78           6.7         3.0          5.0         1.7
101          6.3         3.3          6.0         2.5
103          7.1         3.0          5.9         2.1
104          6.3         2.9          5.6         1.8

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

#or ggplot2
tt = cbind(iris2,a = kmeans.result$cluster)
ggplot(tt,aes(x=Sepal.Width,y=Sepal.Length,colour=a))+geom_point()

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

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

在上图中，聚类中心被标记为星号，异常值标记为’+’

对时间序列进行异常检测

本部分讲述一个对时间序列数据进行异常检测的例子。在本例中，时间序列数据首次使用stl()进行稳健回归分解，然后识别异常值。STL的介绍，请访问 http://cs.wellesley.edu/~cs315/Papers/stl%20statistical%20model.pdf.

#use robust fitting
> f <- stl(AirPassengers,'periodic',robust=TRUE)
> (outliers <- which(f$weights<1e-8))
 [1]  79  91  92 102 103 104 114 115 116 126 127 128 138 139
[15] 140
> 

#set layout
op <- par(mar=c(0,4,0,3),oma = c(5,0,4,0),mfcol=c(4,1))
plot(f,set.pars = NULL)

sts <- f$time.series
#plot outliers
points(time(sts)[outliers],0.8*sts[,"remainder"][outliers],pch = "x",col="red")
par(op)#reset layout

在上图中，异常值用红色标记为’x’

讨论

LOF算法擅长检测局部异常值，但是它只对数值数据有效。Rlof包依赖multicore包，在Windows环境下失效。对于分类数据的一个快速稳定的异常检测的策略是AVF(Attribute Value Frequency)算法。

一些用于异常检测的R包包括：

extremevalues包：单变量异常检测

mvoutlier包：基于稳定方法的多元变量异常检测

outliers包：对异常值进行测验