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,2
1
,+-+---=----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","Test Stat.", "Critical Val.") #给newdf中的列命名
newdf #输出newdf。