R：运行标准差答案

【问题标题】：R: running standard deviationsR：运行标准差
【发布时间】：2012-09-11 13:24:03
【问题描述】：

我需要计算均值向量的“运行”标准差，就像样本量在增加一样计算它们。也就是说，我需要计算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 中可能存在某种基本错误。谢谢！！斯蒂法诺

【问题讨论】：

我会选择se <- sapply(1:length(y),function(x) sd(y[1:x])) 或类似的东西
我认为您经常通过i 潜水。您在se 中除以sqrt(i)，在error 公式中再除以

标签： r montecarlo

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

希望这一切都会有所帮助！

【讨论】：