我想跟进一个关于包中的all.moments
函数的问题moments
:all.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() 数据,那么我可以将其读回并进行操作,而不会消失在我身上。