计算R中for循环的均值

Calculating the means in a for loop in R

大家下午好,

我试图找到我通过模拟 运行 生成的时间序列的标准预测误差,该模拟是通过名为 sim_11 的函数定义的,具有 250 次模拟。这在下面的第一批代码中提供。

第二批创建时间序列模型 (AR(1)) 并尝试预测接下来的 5 个值,我总共进行了 250 次模拟。对于每个模拟,我应该能够得到 5 个预测错误,并且在 250 次模拟之后我应该得到 250 行和 5 列的结果 table。但是,当我尝试在 for 循环中设置它时,我最终只有 250 个单个值,而实际上我应该以 250 x 5 table/matrix 结束。我认为错误在

pred_error_AR1_100[i]<-table((pre_AR1_100$se[1]),(pre_AR1_100$se[2]),
                           (pre_AR1_100$se[3]),(pre_AR1_100$se[4]),
                           (pre_AR1_100$se[5]), ncol=5) 

部分但是我无法弄清楚格式应该在哪里或应该是什么。

提前致谢。

下面提供了两个代码批次用于复制。

# Setup the simulation run with 100 observations and 250 simulations
sim_11=function(){
  e<-rnorm(200, mean=0, sd=0.2) # Produces 200 white noise values
  Y_t=c(0,0)  # Fills in the first 2 observations as a lag of 2 can be handled
  for (i in 3:length(e)){
    f1<- 0.138+(0.316+0.982*Y_t[i-1])*exp(-3.89*(Y_t[i-1])^2)
    f2<- -0.437-(0.659+1.260*Y_t[i-1])*exp(-3.89*(Y_t[i-1])^2)
    Y_t[i]<-f1*Y_t[i-1]+f2*Y_t[i-2]+e[i]
  }
  Y_t<-Y_t[101:200] # Removes the first 100 observations
  Y_t # Prints the 100 observations
}

lapply(1:250, function(x) sim_11()) # Provides the results of the 250 simulations
x_100_lstar=replicate(250,sim_11()) # Places all results into one matrix
pred_error_AR1_100=0
# controls<-list(gammaInt=c(0.1,2000), nGamma=50)
for (i in 1:ncol(x_100_lstar)){
  AR1_100<-ar(x_100_lstar[,i])
  pre_AR1_100<-predict(AR1_100, n.ahead=5)
  pred_error_AR1_100[i]<-table((pre_AR1_100$se[1]),(pre_AR1_100$se[2]),
                           (pre_AR1_100$se[3]),(pre_AR1_100$se[4]),
                           (pre_AR1_100$se[5]), ncol=5)
}
pred_error_AR1_100

要使循环正常工作,您需要将 pred_error_AR1_100 初始化为 n×5 矩阵,然后一次修改一行。你不应该在这里使用 table 。有关构造、访问和修改矩阵的详细信息,请参阅 ?matrix?Extract

n <- ncol(x_100_lstar)
pred_error_AR1_100 <- matrix(NA, n, 5)
for (i in seq_len(n)) {
  AR1_100 <- ar(x_100_lstar[, i])
  pre_AR1_100 <- predict(AR1_100, n.ahead = 5)
  pred_error_AR1_100[i, ] <- pre_AR1_100$se
}

不过,在这些情况下,使用 apply 比自己编写循环更安全、更快速:

## Here, 'x' represents the result of one realization of 'sim_11()'
f <- function(x) {
  AR1_100 <- ar(x)
  pre_AR1_100 <- predict(AR1_100, n.ahead = 5)
  pre_AR1_100$se
}

## Apply function 'f' to each column of 'x_100_lstar'
pred_error_AR1_100 <- t(apply(x_100_lstar, 2, f))

在最后一行,apply的结果,一个5-by-n矩阵,被转置得到一个n-by-5矩阵.

FWIW,如果将 Y_t 初始化为长度为 200 的向量,sim_11() 会稍微快一些,如下所示:

Y_t <- rep.int(NA, 200)
Y_t[1:2] <- 0

而不是在每次迭代中将长度递增 1。