在 R 的 for 循环中计算均值

Posted 2023-02-14

【中文标题】在 R 的 for 循环中计算均值

【英文标题】：Calculating the means in a for loop in R

【发布时间】：2022-01-05 19:46:47

【问题描述】：

大家下午好，

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

第二批创建时间序列模型 (AR(1)) 并尝试预测接下来的 5 个值，我总共进行了 250 次模拟。对于每次模拟，我应该能够得到 5 个预测错误，并且在 250 次模拟之后，我应该有一个 250 行和 5 列的结果表。但是，当我尝试在 for 循环中设置它时，我最终只得到了 250 个单个值，而实际上我应该得到一个 250 x 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) 

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

提前谢谢你。

下面提供了两个代码批处理以供复制。

# 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

【参考方案1】：

要使循环正常工作，您需要将pred_error_AR1_100 初始化为n-by-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 矩阵。

如果您将Y_t 初始化为长度为 200 的向量，那么 FWIW，sim_11() 会稍微快一些，如下所示：

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

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

【讨论】：

嘿 Mikael，谢谢你，因为我之前没有考虑过构造矩阵然后用值填充它。在这段代码的另一段，我做了类似的事情，但是我一次只提取一个值，我可以直接在循环内创建一个 250 x 1 的矩阵。但我知道在 for 循环之外预先构建矩阵而不是在 for 循环内创建它是一个好主意。再次感谢您的帮助。

一般来说，使用apply 比遍历矩阵的列更安全/更快。我已经更新了代码。