在 R 中用循环和 if else 解决 SIR 模型
in R desolve SIR model with loop and if else
我有一个简单的 SIR 模型,我正在尝试实施疫苗接种方法 (V),首先检查感染者是否高于阈值 (100),以及是否仍有足够的易感者 ( 50) 它会在每个时间步 (50) 接种一定数量的疫苗。
但是我想做的是,一旦条件满足,它应该接种 7 天(无论在这 7 天内感染者是否仍然高于阈值,例如,如果在一天之后4、I = 70还是应该继续,只有S < 50才应该停止。7天结束后,应该再次检查条件,或者再开始7天与否。
到目前为止,我所拥有的,
如果有人帮助我实现该循环,我将不胜感激
sirV=function(time, y, params){
S = y[1]
I = y[2]
R = y[3]
V = y[4]
with(as.list(params),{
vac_helper = if (I > 100 & S > 50) {50}
else {0}
N = S+I+R+V
dS = -S*beta*I/N - vac_helper
dI = S*beta*I/N - gamma*I
dR = +gamma*I
dV = vac_helper
return(list(c(dS, dI, dR, dV)))
})
}
myparameters = c(gamma=1/10,beta=0.2)
times <- seq(0, 300)
my_ode <- as.data.frame(ode( y=c(100000, 10, 0,0), times, sirV, myparameters))
这是一个建议,但我不完全确定它是否符合您的要求,所以请先检查一下,如果不符合请再回来找我。
请注意,end
必须在全局环境中。
library(deSolve)
sirV = function(time, y, params){
S = y[1]
I = y[2]
R = y[3]
V = y[4]
with(as.list(params),{
# Has the previous vaccination ended, and do we still need to vaccinate?
if(time > end & I > 100) {
end <<- time + 7
}
# Can we vaccinate?
vaccinate = ifelse(end > time & S > 50, 50, 0)
N = S+I+R+V
dS = -S*beta*I/N - vaccinate
dI = S*beta*I/N - gamma*I
dR = gamma*I
dV = vaccinate
# Store some results in a global list so we can check what is happening under the hood
temp = data.frame(time, end, S, I, R, V, dS, dI, dR, dV, vaccinate)
catch_results[[length(catch_results)+1]] <<- temp
return(list(c(dS, dI, dR, dV)))
})
}
myparameters = c(gamma = 1/10, beta = 0.2)
times <- seq(from = 0, to = 300, by = 1)
end <- 0 # Reset end time
catch_results = list() # Catch results from inside the function
my_ode <- ode( y=c(100000, 10, 0, 0), times, sirV, myparameters)
plot(my_ode)
# Check this to see if we get the expected behavior, especially at around time = 189
results = dplyr::bind_rows(catch_results)
由 reprex package (v0.3.0)
于 2021-08-22 创建
我有一个简单的 SIR 模型,我正在尝试实施疫苗接种方法 (V),首先检查感染者是否高于阈值 (100),以及是否仍有足够的易感者 ( 50) 它会在每个时间步 (50) 接种一定数量的疫苗。
但是我想做的是,一旦条件满足,它应该接种 7 天(无论在这 7 天内感染者是否仍然高于阈值,例如,如果在一天之后4、I = 70还是应该继续,只有S < 50才应该停止。7天结束后,应该再次检查条件,或者再开始7天与否。
到目前为止,我所拥有的, 如果有人帮助我实现该循环,我将不胜感激
sirV=function(time, y, params){
S = y[1]
I = y[2]
R = y[3]
V = y[4]
with(as.list(params),{
vac_helper = if (I > 100 & S > 50) {50}
else {0}
N = S+I+R+V
dS = -S*beta*I/N - vac_helper
dI = S*beta*I/N - gamma*I
dR = +gamma*I
dV = vac_helper
return(list(c(dS, dI, dR, dV)))
})
}
myparameters = c(gamma=1/10,beta=0.2)
times <- seq(0, 300)
my_ode <- as.data.frame(ode( y=c(100000, 10, 0,0), times, sirV, myparameters))
这是一个建议,但我不完全确定它是否符合您的要求,所以请先检查一下,如果不符合请再回来找我。
请注意,end
必须在全局环境中。
library(deSolve)
sirV = function(time, y, params){
S = y[1]
I = y[2]
R = y[3]
V = y[4]
with(as.list(params),{
# Has the previous vaccination ended, and do we still need to vaccinate?
if(time > end & I > 100) {
end <<- time + 7
}
# Can we vaccinate?
vaccinate = ifelse(end > time & S > 50, 50, 0)
N = S+I+R+V
dS = -S*beta*I/N - vaccinate
dI = S*beta*I/N - gamma*I
dR = gamma*I
dV = vaccinate
# Store some results in a global list so we can check what is happening under the hood
temp = data.frame(time, end, S, I, R, V, dS, dI, dR, dV, vaccinate)
catch_results[[length(catch_results)+1]] <<- temp
return(list(c(dS, dI, dR, dV)))
})
}
myparameters = c(gamma = 1/10, beta = 0.2)
times <- seq(from = 0, to = 300, by = 1)
end <- 0 # Reset end time
catch_results = list() # Catch results from inside the function
my_ode <- ode( y=c(100000, 10, 0, 0), times, sirV, myparameters)
plot(my_ode)
# Check this to see if we get the expected behavior, especially at around time = 189
results = dplyr::bind_rows(catch_results)
由 reprex package (v0.3.0)
于 2021-08-22 创建