1

サンプルサイズが増加しているかのように、平均ベクトルの「実行中の」標準偏差を計算する必要があります。つまり、mean_y[1]、mean_y[1] と mean_y[2]、mean_y[1]、mean_y[2]、mean_y[3] などの sd を計算する必要があります。

どこ：

y <- rnorm(10000, mean=5, sd=2) 
i <- seq(1 : length(y))

mean_y <- cumsum(y)/i

同じ基準を適用しようとして (したがって、明示的にループを使用しません)、実行平均のベクトルの標準偏差に対して次のコードを作成しました。

se <- sqrt(1/i^2*cumsum(y^2) - 1/i*mean_y^2)

これは、var(mean(x)) = 1/n*var(x) であるためです。私には、コードは問題ないようです。しかし、i (増加するサンプル サイズ) に対する信頼区間を使用して移動平均をプロットすると、95% のバンドが平均と完全に一致します。

コードは次のとおりです。

error <- qnorm(0.975)*se/sqrt(i)
lower <- mean_y - error
upper <- mean_y + error

# plotting means and ci's against sample size (= up to 10000)

plot(x=i, y=mean_y, xlab="Number of iterations (sample size)", 
ylab="E[y] estimates and 95% CI's", cex=0.4, ylim=c(4.6, 5.4))
lines(lower, col="gold")
lines(upper, col="gold")

理論的根拠は、サンプル サイズが増加し続けるにつれて、推定量 "mean_y" の収束を示すグラフを作成することです。

誰でも私を助けることができますか？seおそらく、式またはlowerとのいずれかに、ある種の基本的なエラーがありupperます。ありがとう！！ステファノ

4

1 に答える 1

1

私が誤ってiで割っていたことが判明しました。また、分散の「n-1」補正を考慮していないため、私の式は Carl Witthoft によって提案されたものと比較できませんでした。

次の行では、最初の問題を解決し、同じグラフをプロットするための 3 つの同等の方法を見つけます。

i <- seq(1 : length(y))
m <- cumsum(y)/i

se_y <- sqrt((1/(i-1)*cumsum(y^2) - i/(i-1)*m^2))

error <- qnorm(0.975)*se_y/sqrt(i)
lower <- m - error
upper <- m + error

# equivalent (slightly slower) methods for getting the std. errors

# method2:
se_2 <- rep(NA, length(y))
for (n in 1:length(y))  {
  se_2[n] <- sd(y[1:n])
}
# method3:
se_3 <- sapply(1:length(y), FUN= function(x) sd(y[1:x]))

最終的なグラフは次のとおりです。

# plotting means and ci's against sample size (= up to 10000)
plot(x=i, y=m, xlab="Number of iterations (sample size)", 
title("Convergence of the ENVP's mean"), 
ylab="E[y] estimates and 95% CI's (EUR millions)", cex=0.4, ylim=c(2620, 2665))
lines(lower, col="gold")
lines(upper, col="gold")
legend("bottomright", legend=c("envp's mean", "95% ci"),
cex=0.8, col=c("black", "gold"), lwd=2, lty=1, bty="n")

dev.copy(tiff, file="mc_envp.tiff", height=6, width=6, units="in", res=200)
dev.off(); dev.off()
windows.options(reset=TRUE)

これがすべて役立つことを願っています！

于 2012-09-11T16:43:43.923 に答える