准备工作
假设检验
- EL for one-dimensional mean
- EL for two-dimensional mean
- EL for three-dimensional mean
一维向量均值的经验似然假设检验
1
2
3
| set.seed(3)
x <- rnorm(50)
(el.x <- el.test(x , mu = 0))
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
| ## $`-2LLR`
## [1] 0.2637869
##
## $Pval
## [1] 0.6075303
##
## $lambda
## [1] -0.08258183
##
## $grad
## [1] -3.876559e-14
##
## $hess
## [,1]
## [1,] 38.86236
##
## $wts
## [1] 0.9264078 0.9764125 1.0218379 0.9131210 1.0164338 1.0024939 1.0071041
## [8] 1.1015784 0.9085495 1.1168962 0.9420583 0.9145632 0.9441460 1.0213091
## [15] 1.0127159 0.9752227 0.9270402 0.9491871 1.1124783 1.0167777 0.9544059
## [22] 0.9278014 0.9834541 0.8790277 0.9615318 0.9423302 1.1060060 1.0912008
## [29] 0.9940828 0.9141791 1.0803514 1.0756629 1.0639385 1.0647605 0.9717422
## [36] 1.0618673 1.1203050 1.0031689 0.9251797 1.0701485 1.0694629 0.9750023
## [43] 1.1631924 0.9384217 1.0296271 0.8424025 0.9867819 1.1030088 0.9637441
## [50] 0.9308779
##
## $nits
## [1] 4
|
此时的x是一维向量, mu 也是长度为 1 的常数向量, 使用 el.test 函数可以执行假设检验.
el.test 函数的返回结果是一个列表, 其中有几个元素较为重要, 分别是
- -2LLR, -2倍的经验似然比, 如果重复带入同一分布的x, 则-2LLR趋近$\chi^2(1)$
- Pval, P($\chi^2(1)$ > -2LLR), 如果小于阈值, 则拒绝 $H_0: \mu = \mu_0$
- lambda, 拉格朗日乘子
- wts, 每个观测上的权重$wi$, 不过这些权重和服从$\sum{i=1}^n w_i = n$
我们可以简单的对以上返回结果进行验证
以下两项相等
1
| 1 - pchisq(el.x$`-2LLR`, df = 1)
|
以下三项相等
1
| 2*sum(log(1 + el.x$lambda*x))
|
以下两项相等
二维向量均值的经验似然假设检验
对于二维情况, 代码如下:
1
2
3
4
| set.seed(3.1)
x <- rnorm(50)
y <- rnorm(50)
(el.two.dim <- el.test(cbind(x,y), mu = c(0,0)))
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
| ## $`-2LLR`
## [1] 0.7596482
##
## $Pval
## [1] 0.6839817
##
## $lambda
## [1] -0.07377731 0.11984790
##
## $grad
## [1] -2.838477e-10 3.860352e-10
##
## $hess
## x y
## x 39.141566 -3.804552
## y -3.804552 35.404447
##
## $wts
## [1] 0.8634990 1.0815816 0.9872478 1.1405190 1.2248908 1.0599729 1.1499584
## [8] 0.9251811 0.8334172 1.2309218 0.8899308 0.8378703 0.9228783 0.9769991
## [15] 0.8853284 1.0361947 0.9802815 0.8603676 1.3241724 0.9924831 0.9625337
## [22] 0.8884932 0.8789057 0.8658081 0.9402688 0.8739183 0.9431862 1.0295723
## [29] 1.1275207 0.9389052 0.8759984 1.1173978 0.9719654 0.9154111 0.8919223
## [36] 1.0557637 1.0748899 1.1224138 0.8889258 1.1973087 1.1909791 1.1058810
## [43] 1.2552493 0.8439272 1.0782284 0.8630544 1.0453786 1.0189582 0.8731615
## [50] 0.9603777
##
## $nits
## [1] 4
|
将两个向量按列合并为一个矩阵, mu也是一个长度为2的常数向量. 对el.test的返回结果, 我们也有一维情形下的等式
以下两项相等
1
| 1 - pchisq(el.two.dim$`-2LLR`, df = 2)
|
以下三项相等
1
| 2*sum(log(1 + cbind(x,y) %*% el.two.dim$lambda))
|
1
| -2*log(prod(el.two.dim$wts))
|
以下三项相等
三维向量均值的经验似然假设检验
1
2
3
4
5
| set.seed(3.14)
x <- rnorm(100)
y <- rnorm(100)
z <- rnorm(100)
el.three.dim <- el.test(cbind(x,y,z), mu = rep(0,3))
|
将三个向量按列合并为一个矩阵, mu也是一个长度为3的常数向量. 对el.test的返回结果, 我们也有一维情形下的等式
以下两项相等
1
| 1 - pchisq(el.three.dim$`-2LLR`, df = 3)
|
以下三项相等
1
| 2*sum(log(1 + cbind(x,y,z) %*% el.three.dim$lambda))
|
1
| -2*log(prod(el.three.dim$wts))
|
以下四项相等
1
| length(x);length(y);length(z)
|
当 $\mu$ 远离 $\mu_0$ 时, el.test 的 Pval 会变小
1
2
3
4
5
6
7
8
9
10
11
12
| set.seed(3.141)
x <- rnorm(100)
n <- 10
pval <- numeric(n)
minus2LLR <- numeric(n)
mu <- seq(0, by = 0.02, length.out = n)
for(i in 1:n){
el.x <- el.test(x, mu = mu[i])
pval[i] <- el.x$Pval
minus2LLR[i] <- el.x$`-2LLR`
}
data.frame(mu,pval, minus2LLR)
|
1
2
3
4
5
6
7
8
9
10
11
| ## mu pval minus2LLR
## 1 0.00 0.89697380 0.01676640
## 2 0.02 0.91614625 0.01108584
## 3 0.04 0.73346985 0.11594945
## 4 0.06 0.56451728 0.33194288
## 5 0.08 0.41670823 0.65958007
## 6 0.10 0.29441724 1.09931157
## 7 0.12 0.19875088 1.65152955
## 8 0.14 0.12800203 2.31656962
## 9 0.16 0.07854717 3.09470955
## 10 0.18 0.04587538 3.98616465
|
可以看出, 当 $\mu$ 离真值$\mu_0 = 0$ 越来越远时, p值(pval)越来越小, -2LLR 值越来越大
置信区间
经验似然置信区间的求法:
- 先定义一个函数(myfun)
- myfun的第一个参数应该是正在给哪个参数求置信区间
- myfun的返回值应该是一个列表, 并且其中含有一个元素名称为”-2LLR”
- findUL函数可以用来求置信区间的上下限, 通过level参数可以设置置信度
具体使用方法如下:
1
2
3
| myfun <- function(theta, x){
el.test(x, mu = theta)
}
|
1
| findUL(step = 0.1, fun = myfun, MLE = mean(x), level = 3.84, x = x)
|
一维向量均值的经验似然置信区间
1
2
3
| set.seed(3)
x <- rnorm(50)
findUL(step = 0.1, fun = myfun, MLE = mean(x), level = 3.84, x = x)
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| ## $Low
## [1] -0.31201
##
## $Up
## [1] 0.1799538
##
## $FstepL
## [1] 1e-09
##
## $FstepU
## [1] 1e-09
##
## $Lvalue
## [1] 3.84
##
## $Uvalue
## [1] 3.84
|
findUL函数返回结果是个列表, 其中的Low和Up分别是置信区间的下限和上限
需要注意的是, findUL函数目前只能对一维参数进行求置信区间上下限, 而无法对多维参数求置信域
经验似然和正态逼近的对比
经验似然(Empirical likelihood)和正态逼近(Normal Approximation)在对参数进行区间估计时, 一般从置信区间的覆盖率(coverage probability, 简称cp)和区间平均长度(average length, 简称al)来进行比较. 这就需要生成多个置信区间取平均进行比较.
经验似然方法计算置信区间
生成100组服从标准正态分布的随机样本, 每组样本的样本量为50, 可以得到100个置信区间的上下限
1
2
3
4
5
6
7
8
9
10
11
| set.seed(3.1415)
elci.up <- numeric(100)
elci.low <- numeric(100)
elci.length <- numeric(100)
for(i in 1:100){
x <- rnorm(50)
elci <- findUL(step = 0.1, fun = myfun, MLE = mean(x), x = x)
elci.up[i] <- elci$Up
elci.low[i] <- elci$Low
elci.length[i] <- elci.up[i] - elci.low[i]
}
|
我们可以计算这100个置信区间的平均长度(average length of empirical likelihood, 记为alel)
1
| (alel <- mean(elci.length))
|
还可以计算着100个置信区间的覆盖率(coverage probability of empirical likelihood, 记为cpel)
1
| (cpel <- sum(elci.up >= 0 & elci.low <= 0)/100)
|
我们还可以把这100个置信区间作图, 并把没有覆盖到真值的置信区间用红色标注出来
1
2
3
4
5
6
7
8
9
10
| plot(1:100, seq(-1,1,length.out = 100), type = "n",xlab = "x", ylab = "y")
points(1:100, elci.low)
points(1:100, elci.up)
outofci <- which(elci.up < 0 | elci.low > 0)
points(outofci, elci.low[outofci], pch = 16, col = "red")
points(outofci, elci.up[outofci], pch = 16, col = "red")
for(i in 1:100){
segments(i,elci.low[i],i,elci.up[i])
}
segments(0,0,100,0)
|
正态逼近方法计算置信区间
正态逼近使用t检验(t.test)来获取样本的置信区间上下限
1
2
3
4
5
6
7
8
9
10
| set.seed(3.1415)
naci.up <- numeric(100)
naci.low <- numeric(100)
naci.length <- numeric(100)
for(i in 1:100){
x <- rnorm(50)
naci.up[i] <- t.test(x)$conf.int[2]
naci.low[i] <- t.test(x)$conf.int[1]
naci.length[i] <- naci.up[i] - naci.low[i]
}
|
我们可以计算这100个置信区间的平均长度(average length of normal approximation, 记为alna)
1
| (alna <- mean(naci.length))
|
还可以计算着100个置信区间的覆盖率(coverage probability of normal approximation, 记为cpna)
1
| (cpna <- sum(naci.up >= 0 & naci.low <= 0)/100)
|
在目前的比较中, EL和NA的置信区间覆盖率均为0.96. 对于置信区间的平均长度, EL为0.5581938, 而NA为0.5705212.
EL和NA置信区间综合比较
将上述生成置信区间的方法写成函数elna, 并用apply函数替代循环
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
|
elul <- function(x){
findUL(step = 0.1, fun = myfun, MLE = mean(x),
level = 3.84, x = x)
}
elna <- function(repe,size,seed){
set.seed(seed)
time1.length <- Sys.time()
x <- matrix(rnorm(size*repe), nrow = size, ncol = repe)
na <- unlist(apply(x, MARGIN = 2, t.test))
naci.up <- as.numeric(na[seq(5, by = 10, length.out = repe)])
naci.low <- as.numeric(na[seq(4, by = 10, length.out = repe)])
naci.length <- naci.up - naci.low
el <- unlist(apply(x, MARGIN = 2, FUN = elul))
elci.up <- el[seq(2, by = 6, length.out = repe)]
elci.low <- el[seq(1, by = 6, length.out = repe)]
elci.length <- elci.up - elci.low
compare <- data.frame(naci.up, naci.low, naci.length,
elci.up, elci.low, elci.length,
na_el = naci.length - elci.length)
time2.length <- Sys.time()
list(na_el_num = sum(compare["na_el"] > 0),
alna = mean(naci.length),
alel = mean(elci.length),
cpna = sum(naci.up >= 0 & naci.low <= 0)/repe,
cpel = sum(elci.up >= 0 & elci.low <= 0)/repe,
time = time2.length - time1.length)
}
|
如果想构造100个置信区间, 每个置信区间是由样本量为30的样本产生的, 那么设置种子为3, 程序调用如下
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
| ## $na_el_num
## [1] 99
##
## $alna
## [1] 0.7271278
##
## $alel
## [1] 0.699543
##
## $cpna
## [1] 0.96
##
## $cpel
## [1] 0.96
##
## $time
## Time difference of 50.82591 secs
|
可以看到, 程序耗时约为循环次数的三分之一. 在构造的100个置信区间中, 有99个的NACI都比ELCI要长, 而平均长度也能说明EL要优于NA. 从覆盖率角度来看, NACI的覆盖率与ELCI相当.
事实上, 当样本量增大时, 无论是EL还是NA的置信区间长度都在缩短, 而覆盖率也都在提升. 此时两个方法的效果基本相当.
