如何解决运行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}
在第j次仿真中,S_ {i,j} = S [i] =第i天的易感人群。
然后,我可以找到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)
}
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。