奇怪的 all.moments() 结果答案

【问题标题】：Odd all.moments() outcome奇怪的 all.moments() 结果
【发布时间】：2012-12-16 09:27:51
【问题描述】：

我想跟进一个关于 moments 包中的 all.moments 函数的问题：all.moments function weird outcome

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

考虑以下代码：

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() 数据，然后我可以将其读回并进行操作，而不会消失在我身上。

【问题讨论】：

我不会提供证明，但我觉得问题可能在于 all.moments 默认情况下不计算 central 时刻。试试central = TRUE。
我会尝试 - 当随机法线的平均值为 0 时，给定 1M 样本大小，并且考虑到差异有多大，我有点怀疑，但你很可能是对的。跨度>

标签： r statistics

【解决方案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

【讨论】：