Hi,
does the following help you?
Thomas
require(deSolve)
SEIR <- function(t, x, p) {
if (t < 7)
xlag <- x
else
xlag <- lagvalue(t - 7)
V <- ifelse(xlag[3]/sum(xlag) > 0.01, 0.25, 0)
## uncomment the following for printing some internal information
#cat("t=", t, " -- ")
#cat(xlag, " -- ", V, "\n")
N <- sum(x)
with(as.list(c(x,p)),{
dS<- b*N - d*S - beta*S*I/N - V*S
dE<- -d*E+beta*S*I/N - epsilon*E - V*E
dI<- -d*I + epsilon*E - gamma*I - mu*I - V*I
dR<- -d*R + gamma*I + V*S + V*E + V*I
list(c(dS, dE, dI, dR), N = unname(N), V = unname(V))
})
}
num_years <- 1
time_limit <- num_years*365.00
b <- 1/(10.0*365)
d <- b
beta <- 0.48
epsilon <- 1/4
gamma <- 1/4
mu <- -log(1-0.25)*gamma
parms <- c(b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu)
xstart <-c(S=999, E=0, I=1, R=0)
times <- seq(0.0, time_limit, 1.0)
out <- dede(xstart, times, SEIR, parms, rtol=1e-8, atol=1e-8)
plot(out)
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
