0

我想跟进一个关于包中的all.moments函数的问题momentsall.moments 函数奇怪的结果

all.moments 输出的分布 - 特别是渐近方差 - 似乎与第 3 和第 4 时刻的标准理论不一致,或者我做错了什么。

考虑以下代码:

reps <- 100000                              # number of repetitions
n <- 1000000                                # sample size
asy.cmtx <- matrix(0,nrow=reps,ncol=5)      # initialize our matrix of conventional moments
for(i in 1:reps){
vec <- rnorm(n,mean=0,sd=1)               # create vector of 1M standard normals
asy.cmtx[i,] <- all.moments(vec,order.max=5)[2:6]    # get sample moments 
}
mean(asy.cmtx[,3])             # this should be 0
# [1] 3.972406e-06
mean(asy.cmtx[,4])            # this should be 3
# [1] 2.999996
var(sqrt(n)*asy.cmtx[,3]/(asy.cmtx[,2]^(3/2)))     # this should be 6
# [1] 14.41766
var(sqrt(n)*(asy.cmtx[,4]-3)/(asy.cmtx[,2]^(2))) # this should be 24
# [1] 96.25745

渐近方差似乎比预期的要大。

我回去简单地使用了我们都知道的示例矩公式并得到了正确的结果:

reps <- 100000                              # number of repetitions
n <- 1000000                                # sample size
asy.34.2 <- matrix(0,nrow=reps,ncol=2)      # initialize our matrix of conventional moments
for(i in 1:reps){
y <- rnorm(n,mean=0,sd=1)               # create vector of 1M standard normals
y.bar <- mean(y)
m2 <- ((1/n)*sum((y-y.bar)^2))
asy.34.2[i,1] <- ((1/n)*sum((y-y.bar)^3))/(m2^(3/2))
asy.34.2[i,2] <- ((1/n)*sum((y-y.bar)^4))/(m2^2)
}

mean(asy.34.2[,1])
# [1] 7.512593e-06
mean(asy.34.2[,2])
# [1] 2.999985
var(sqrt(1000000)*asy.34.2[,1])     # this should be 6
# [1] 5.990771
var(sqrt(1000000)*(asy.34.2[,2]-3)) # this should be 24
# [1] 24.23367                    

所以实际的公式有我们期望的方差,但 all.moments 没有,即使它有正确的期望。由于样本大小或迭代次数,这应该不是问题,因为在这种情况下两者都很大(两种情况下分别为 1M 和 100k)。以前有人遇到过这个问题吗?

另外,有没有人对 R 简单地丢弃你的输出有问题?当我试图操纵上面的输出时,输出以某种方式“消失”了,dim = NULL。但是,如果我做的第一件事是 write.csv() 数据,那么我可以将其读回并进行操作,而不会消失在我身上。

4

1 回答 1

1

有两个错误:(1)正如弗洛德尔指出的那样,您需要中心矩而不是原始矩,这对矩的点估计没有太大区别,但对它们的方差估计确实有影响;(2) 使用 all.moments 计算峰度时存在错误。当我们应用这些更正时,我们有:

reps <- 1000                              
n <- 100000

...

    asy.cmtx[i,] <- all.moments(vec,order.max=5, central = TRUE)[2:6]

...

var(sqrt(n)*asy.cmtx[,3]/(asy.cmtx[,2]^(3/2)))     # this should be 6
# [1] 5.993772
var(sqrt(n)*(asy.cmtx[,4])/(asy.cmtx[,2]^(2)) - 3) # this should be 24
# [1] 23.89204
于 2013-01-02T06:29:31.757 回答