使用for循环保存模拟结果

Question

我正在尝试运行更复杂的模拟，并将reach复制的结果保存到数据帧中。我想要做的一个非常简单的抽象是这样的：

sim.res <- data.frame(mn1 = mn1, mn2 = mn2, rep = rep)

for (k in 1:10){
  
  set.seed(k*8)
  mn1 <- mean(rnorm(25, 1, 1))
  mn2 <- mean(rnorm(25, 2, 1))

  sim.res$mn1 <- mn1
  sim.res$mn2 <- mn2
  sim.res$rep <- k
}

输出应该是一个具有三列的数据帧，在本例中为10行，每个复制的结果存储在每行中。我得到的是第10次复制。我哪里错了？

谢谢。

Answer 1

每次迭代循环时都会覆盖每一列的完整性，因为默认情况下R是矢量化的(在列上运行)。而且，没有必要在每次迭代循环时重置种子。试试这样的东西。

library(tidyverse)

set.seed(123)

lapply(
  1:10,
  function(k) {
    data.frame(
      mn1 = mean(rnorm(25, 1, 1)), 
      mn2 = mean(rnorm(25, 2, 1)), 
      rep = k
    )
  }
) %>%
bind_rows()

         mn1      mn2 rep
1  0.9666697 2.102137   1
2  1.0102410 2.282576   2
3  0.7231910 1.769008   3
4  1.2152427 1.862371   4
5  1.1074981 1.875461   5
6  1.1727220 2.326180   6
7  0.9789770 2.025131   7
8  1.0981718 1.752828   8
9  0.9402713 1.928519   9
10 1.2714378 2.283175  10

或者，利用从N(µ，sd)得出的n的平均值是N(mu，sd/sqrt(n))的事实，

data.frame(
  rep=rep(1:10),
  mn1=rnorm(10, 1, sqrt(1/25)),
  mn2=rnorm(10, 2, sqrt(1/25))
)

   rep       mn1      mn2
1    1 1.2053570 2.176493
2    2 1.1502123 2.041120
3    3 0.6981667 1.876713
4    4 0.9809705 1.853040
5    5 0.8208104 1.973639
6    6 0.5858498 2.062003
7    7 1.0300240 1.792064
8    8 0.9841577 1.963138
9    9 0.9805261 2.193453
10  10 1.0432305 1.978344

Answer 2

在R中，最基本的对象仍然是一个向量，所以这根本不需要循环。

当您设置sim.res$mn1 <- mn1时，您将把mn1放入data.frame中，但是由于mn1的长度仅为1，而sim.res有10行，因此R将该值重复10次并将其放入data.frame中。

此外，为了重现性，set.seed只需要使用一次。

set.seed(42)
sim.res <- data.frame(mn1 = replicate(10, mean(rnorm(25, 1, 1))),
                      mn2 = replicate(10, mean(rnorm(25, 2, 1))), 
                      rep = 1:10)
sim.res
         mn1      mn2 rep
1  1.0818928 1.934759   1
2  1.0153544 2.185476   2
3  0.9898056 1.523066   3
4  0.9132660 2.229027   4
5  1.1003520 1.754576   5
6  0.7468452 2.211500   6
7  0.8480360 2.303661   7
8  1.2965247 2.014611   8
9  0.7511059 1.801408   9
10 1.0231177 1.924155  10

replicate只是sapply的一个包装器，它使得这段代码更容易理解。

也许从技术上讲，这是在使用循环，但它看起来肯定不像循环。

Answer 3

感谢您的回复。我还找到了一个解决方案，它给出了我想要的结果：

sim.res <- list() # List container

for (k in 1:10){
  mn1 <- mean(rnorm(25, 1, 1))
  mn2 <- mean(rnorm(25, 2, 1))
  
  tmp <- cbind(mn1, mn2, k)
  
  sim.res[[k]] <- tmp
  
}

(我知道这不是最有效的)

最好的，莱斯利