ESD—检验离群值

Generalized ESD Test for Outliers 一、广义ESD 检验是做什么的

广义ESD 检验是一个检测离群值的方法。它检验服从近似正态 分布的一个单变量数据集中的一个或多个离群值。 二、为什么要使用这个算法

许多统计技术对离群值的存在是敏感的。例如，计算一个数据集 的均值或标准差时，离群值的影响是很大的。因此，检验离群值应该

是任何数据分析的常规部分。我们对潜在的异常值进行检查，以查看

它们是否可能是错误的。如果数据点是错误的，但如果可能，应当校

正，如果不可能则删除。如果没有理由相信边远点是错误的，它不应

该在没有仔细考虑的情况下被删除。 三、对广义ESD 检验的定义

给定数据集X=（x1，x2,...,xn ），设),(~2 σμN X ，x1，

x2,...,xn 相互独立且与X 有相同的概率分布。首先画出数据 集的正态概率图（运行序列图，箱线图，或直方图），观察是否存在

潜在离群值（若事先不知道数据是否服从近似正态分布，还可评估数

据是否遵循一个近似正态分布），以确定是否有必要进行离群值检验。

若存在离群值，则给定的离群值数目的上限，令为r ，则广义ESD 检

验实质上是执行r 次单独的检验：首先检验第一个可能的离群值，

计

算相应的统计量，在给定的显著水品α下做出判断；再检验第二个离群值，...，检验第r个离群值。这r次检验相互独立，互不影响。具体地说，我们假设：

H0 ：没有离群值 H1 ：最多有r个离群值 计算检验统计量Ri ： s x x Ri i | | max- =

x：表示样本均值s：表示样本标准差 公式中 | |x

x i-的值越大，说明i x与x相差越大，该数距点是 离群值的可能性也越大。我们首先删除使 | |x

x i-最大的i x，然

后重新计算余下的n-1个数据的Ri，再移除相应的i

x，重复这个过程，一直到移除了r个满足条件的数据（此时，该数据集中，可能是离群值的r个数据被删除），形成r检验统计量R1，R2，...,Rr。

在显著性水平为α（置信度为1-α）的条件下，计算检验的临界值i

λ

)1)(1()(1,21

,+-+---=----i n i n t i n i i n p i n p t λ

1,--i n p t ：自由度为n-i-1的t 分布的100p 百分位点 而)1(21+--=i n p α

由相关知识知道t 的密度函数 为 )11(2)2

1()1()2()()2(--+----Γ?---Γ=i n t i n i n i n i n t f π 假设H0成立，则有?=i p dt t f λ0 )( 则有 p i t P =≤)(λ

则上述检验的拒绝域为),(+∞i λ，即当Ri >i λ时，对应的数据是离群值。从而找出使i R > i λ得最大的i ，就是我们检验的数据集中存在i 个离群值。

注意：研究表明，当n≧25时，临界值是非常正确的；当n≧15时，临界值是较正确的。

四、具体算法

y <- c(-0.25, 0.68, 0.94, 1.15, 1.20, 1.26, 1.26, 1.34, 1.38, 1.43, 1.49, 1.49, 1.55, 1.56, 1.58, 1.65, 1.69, 1.70, 1.76, 1.77, 1.81, 1.91, 1.94, 1.96, 1.99, 2.06, 2.09, 2.10,

2.14, 2.15, 2.23, 2.24, 2.26, 2.35, 2.37, 2.40, 2.47, 2.54, 2.62, 2.64, 2.90, 2.92, 2.92, 2.93,

3.21, 3.26, 3.30, 3.59, 3.68, 4.30, 4.64, 5.34, 5.42,

6.01) #输入待监测的数据

qqnorm(y) #画出y的正态概率图，观察是否存在潜在离群值 评估数据是否遵循一个近似正态分布 #构建函数计算检验统计量r rval = function(y){

ares = abs(y - mean(y))/sd(y) #计算y与y均值差的绝对值与y

的标准差的商

df = data.frame(y, ares) #创建数据框，包含两列：y ares r = max(df$ares) #找出ares的最大值并赋给r

list(r, df) #建一个包含r与df的列表} ## 定义数值和向量

n = length(y) #将y的长度值赋予n alpha = 0.05 #显著性水平为0.05 lam = c(1:10) #为lam赋初值 R = c(1:10) #为R赋初值

## 计算检验统计量，每次循环，删除当前数据集中与该数据集的均值相差最大的数值

for (i in 1:10){ if(i==1){ #第一次循环

rt = rval(y) #对y调用函数rval，将其返回值赋给rt

R[i] = unlist(rt[1]) #将rt的第一个对象r（即上述ares的最大值）从列表中取出来，并将其类型转变为以一个普通的矢量,接着赋给R[i]

df = data.frame(rt[2])#取出rt中的第二个对象，将其封装到数据框df中

newdf = df[df$ares!=max(df$ares),]}# 将df中zres最大的行删除，其余的行赋给newdf

else if(i!=1){#除第一次循环外都执行

rt = rval(newdf$y)#对newdf中的对象y调用函数rval，将返回值赋予rt

R[i] = unlist(rt[1]) df = data.frame(rt[2])

newdf = df[df$ares!=max(df$ares),]} ## 计算临界值lam p = 1 - alpha/(2*(n-i+1))

t = qt(p,(n-i-1)) #算出自由度为（n-i-1）的t分布的100p分为点

lam[i] = t*(n-i) / sqrt((n-i-1+t**2)*(n-i+1)) }

## 输出结果

newdf = data.frame(c(1:10),R,lam) #将向量（1:10），R，lam封装到数据框newdf中

names(newdf)=c(\"No. Outliers\给newdf中的列命名

newdf #输出newdf