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公式中还是在lowerand中upper。谢谢！！斯特凡诺

4

1 回答 1

1

事实证明，我错误地除以 i。此外，由于不考虑方差的“n-1”校正，我的公式与 Carl Witthoft 建议的公式不具有可比性。

在以下几行中，您可以找到解决初始问题并绘制相同图形的三种等效方法：

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 回答