r - 运行模拟 n 次并将值存储在矩阵中并取平均值
问题描述
我正在对流行病数据进行模拟,我计算了 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
解决方案
如果您想将数据保存在矩阵中以进行多次模拟,这里是一个示例
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)
}
推荐阅读
- android - 制作带有菜单图标的工具栏,如亚马逊应用程序
- wpf - 突出显示 DataGrid 中的空单元格
- php - 拉拉维尔 5.6。如何使用 html 锚点或其他标签切换语言?
- css - Vuetify - v-btn(按钮)的标签与横向边界重叠
- kentico - 如何找到跨 Kentico 实例使用表单的每个页面?
- java - 工具栏 findViewById 返回空指针
- node.js - 如何终止 npm 查询器提示并将控制权返回到主菜单/功能
- linux - GNU makefile shell 命令无法正常工作
- mysql - 如何检查MySQL中是否存在主键
- python - Python 请求:登录 Microsoft.com 但退出