运行模拟 n 次并将值存储在矩阵中并取平均值

Question

我正在对大流行数据进行模拟，我计算了 MLE，其值例如为 0.99。这是用于 SEIR 建模的，所以我有一个 S、E、I 和 R 的数据框。现在我正在运行相同的模拟，但我想复制模拟 100 次，然后考虑平均值。

我的模拟代码，如下。

### Pre-define VALUES
# 50 days
sumofnew_infec<-rep(0,50)
Snew<-rep(0,50)
Enew<-rep(0,50)
Inew<-rep(0,50)
Rnew<-rep(0,50)
Snew[1]<-Current_dayStats$St[1]
Inew[1]<-Current_dayStats$It[1]
Enew[1]<-Current_dayStats$Et[1]
Rnew[1]<-Current_dayStats$Rt[1]

E_I<-0
I_R<-0

### SIMULATION STARTS HERE
for(i in 1:49)
{
  newinfections<-rbinom(n=Snew[i],size=1,prob=(1-MLE^Inew[i]))
  sumofnew_infec[i]<-sum(newinfections)
  Snew[i+1]<-Snew[i]-sumofnew_infec[i]
    if(i>0)
    {  
      E_I<-sum(sumofnew_infec[i])
      #E_I<-0
      I_R<-sum(Enew[i])
    }
    else
    {
      E_I<-sumofnew_infec[i]
      I_R<-sum(Enew[i])
    }
    Enew[i+1]<-Enew[i]-sumofnew_infec[i]+E_I
    Inew[i+1]<-Inew[i]+E_I-I_R
    Rnew[i+1]<-I_R+Rnew[i]
}
sumofnew_infec
Snew
Enew
Inew
Rnew

我想将结果存储在一个矩阵中，例如

S = S_{i,j}

其中 S_{i,j} = S[i] = 第 i 天的易感人群，在第 j 次模拟中。

然后我可以找到 S_{i,1}, S_{i,2}, ..., S_{i,100} 的平均值，这将是第 i 天易感者数量的平均模型预测。最后，我可以绘制所有这些平均值以查看平均易感过程。这就是整体，我正在尝试使用复制，在函数之上创建，但这不起作用。任何帮助，将不胜感激。提前致谢。

编辑： 我在一个函数中创建了模拟。

> do_once()
 [1] 180 176 173 167 155 136 105  57  19   3   1   0   0   0   0   0   0   0   0   0   0   0   0
[24]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
[47]   0   0   0   0

Answer 1

如果您想将数据保存在矩阵中以进行多次模拟，这里是一个示例

dayNum <- 50
simNum <- 100
S <- matrix(dayNum*simNum,nrow = dayNum)
for (j in 1:simNum) {
  for (i in 1:dayNum) {
    S[i,j] <- runif(1)
  }
}

---

当你想计算S的平均值时，你可以使用rowMeans

rowMeans(S)

---

编辑

sumofnew_infec_out <- c()
Snew_out <- c()
Enew_out <- c()
Inew_out <- c()
Rnew_out <- c()
for (k in 1:100) {
  for (i in 1:49)
  {
    newinfections <- rbinom(n = Snew[i], size = 1, prob = (1 - MLE^Inew[i]))
    sumofnew_infec[i] <- sum(newinfections)
    Snew[i + 1] <- Snew[i] - sumofnew_infec[i]
    if (i > 0) {
      E_I <- sum(sumofnew_infec[i])
      # E_I<-0
      I_R <- sum(Enew[i])
    }
    else {
      E_I <- sumofnew_infec[i]
      I_R <- sum(Enew[i])
    }
    Enew[i + 1] <- Enew[i] - sumofnew_infec[i] + E_I
    Inew[i + 1] <- Inew[i] + E_I - I_R
    Rnew[i + 1] <- I_R + Rnew[i]
  }
sumofnew_infec_out <- cbind(sumofnew_infec_out,sumofnew_infec)
Snew_out <- cbind(Snew_out,Snew)
Enew_out <- cbind(Enew_out,Enew)
Inew_out <- cbind(Inew_out,Inew)
Rnew_out <- cbind(Rnew_out,Rnew)
}