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

2つのエラーがあります。(1)flodelが、生のモーメントではなく中心モーメントが必要であると指摘しているため、これはモーメントの点推定には大きな違いはありませんが、分散の推定には違いがあります。(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 に答える