r - 未将行添加到R中的数据框
问题描述
所以,我一直在尝试对 SIR 模型进行模拟研究。
我有以下代码(试图清理它):
# Initial parameters
N <- 1E6 # Total population
I <- 1 # Number of Infectious at time 0
S <- N-1 # Number of Susceptibles at time 0
R <- 0 # Number of Recovered at time 0
# Vector to store observations in
Df1 <- data.frame("final_size" = as.numeric(), "peak_size" = as.numeric())
# Setting a seed for reproducibility
set.seed(1996)
n_sim <- 100
# Setting different values for R0, the basic reproduction number
R0 <- seq(0.5,2.5, 1)
for(values in R0){
if(values == 0.5){
# Transmission parameters
R0 <- values # Basic Reproduction number
nu <- 1/6 # Recovery rate (in days)
b <- nu*R0/N # Infection rate (in days)
for(sim in n_sim){
temp <- NULL
# Binomial model
#----------------
# Initial states
Sold = S # Number of Susceptibles at time t=0
Iold = I # Number of Infectious at time t=0
Rold = R # Number of Recovered at time t=0
# Output vectors
Svec =Sold; Ivec = Iold; Rvec = Rold
stop = FALSE
# Loop - continue until stop=TRUE
while (!stop){
Ih = rbinom(1,Sold,(1-exp(-b*Iold)))
Rh = rbinom(1,Iold,(1-exp(-nu)))
Rh = nu*Iold
Sold = Sold-Ih
Iold = Iold+Ih-Rh
Rold = Rold+Rh
Svec = c(Svec,Sold)
Ivec = c(Ivec,Iold)
Rvec = c(Rvec,Rold)
if (Iold<=2e-5){stop=T}
}
peak_size_df <- max(Ivec)
final_size_df <- Rvec[length(Rvec)]/N
temp <- rbind(temp, c(final_size_df, peak_size_df))
colnames(temp) <- c("final_size", "peak_size")
Df1 <- rbind(Df1, temp)
}
}
}
我正在寻找将数据存储在 Df1 中。但是,在循环结束时,只存储了 1 个循环,我假设是最后一个循环。我不太明白。我已经重新编码了几次,在这些情况下,我最终得到了 99 NA 行的最终大小和峰值大小。在这个版本中,我只得到 1 行(但是有值)。我计划为 R0 序列中看到的不同 R0 值扩展循环。由于它不适用于第一个值,我还没有扩展它。有什么建议么?改进?
在 Gregory 和 Cath 发表意见后,进行以下调整:
# Setting a seed for reproducibility
#set.seed(1996)
n_sim <- 100
# Initial parameters
N <- 1E6 # Total population
I <- 1 # Number of Infectious at time 0
S <- N-1 # Number of Susceptibles at time 0
R <- 0 # Number of Recovered at time 0
# Vector to store observations in
Df1 <- data.frame("final_size" = rep(NA, n_sim), "peak_size" = rep(NA, n_sim))
Df2 <- NULL
Df3 <- NULL
# Setting different values for R0, the basic reproduction number
R0 <- seq(0.5,2.5, 1)
#plot(Svec, type = "l", ylim = c(0, 1000000), col = "red")
#lines(Rvec, type = "l", col = "blue")
#lines(Ivec, type = "l")
#max(Ivec)
for(values in R0){
if(values == 0.5){
# Transmission parameters
R0_value <- values # Basic Reproduction number
nu <- 1/6 # Recovery rate (in days)
b <- nu*R0_value/N # Infection rate (in days)
for(sim in n_sim){
# Binomial model
#----------------
# Initial states
Sold = S # Number of Susceptibles at time t=0
Iold = I # Number of Infectious at time t=0
Rold = R # Number of Recovered at time t=0
# Output vectors
Svec =Sold; Ivec = Iold; Rvec = Rold
stop = FALSE
# Loop - continue until stop=TRUE
while (!stop){
Ih = rbinom(1,Sold,(1-exp(-b*Iold)))
Rh = rbinom(1,Iold,(1-exp(-nu)))
Rh = nu*Iold
Sold = Sold-Ih
Iold = Iold+Ih-Rh
Rold = Rold+Rh
Svec = c(Svec,Sold)
Ivec = c(Ivec,Iold)
Rvec = c(Rvec,Rold)
if (Iold<=2e-5){stop=T}
}
peak_size_df <- max(Ivec)
final_size_df <- Rvec[length(Rvec)]/N
Df1[sim, "final_size"] <- final_size_df
Df1[sim, "peak_size"] <- peak_size_df
}
}
}
这给了我 99 行 NA 和第 100 行最后一个模拟值。知道是什么导致了NA吗?
解决方案
推荐阅读
- scala - spark知道在kafka中添加的新分区吗?
- sql - PostgreSQL:在最近一次考试之日参加过一项以上考试的学生
- multipartform-data - 如何在 Spring REST 中将 json 响应拆分为具有给定边界值的 MultipartForm 数据
- php - 使用 Cybersource SOAP API 时如何获取付款令牌
- c# - 无法等待异步任务
- symfony - Symfony 4 twig 学说,将子菜单关联到下拉列表
- sql - SQL 查询以在一行中获取节点的所有祖先
- sharepoint - 搜索 Refiner & Refinable 字符串
- python - keras LSTM 预测不佳
- c# - 此解析操作已结束。Autofac、Automapper 和 IMemberValueResolver