Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup--------------------------------------------------------------------
library(MicroMoB)
library(ggplot2)
library(data.table)
library(parallel)
## ---- echo=FALSE--------------------------------------------------------------
pardefault <- par()$mar
## ---- fig.width=8, fig.height=6-----------------------------------------------
tmax <- 365 * 3
p <- 3
lambda <- dnorm(x = 1:365, mean = 180, sd = 90)
lambda <- lambda * (100/max(lambda))
lambda <- t(replicate(p, lambda))
psurv <- (sin((1:365)/365*2*pi) + 1.01)/2 * 0.9
EIP <- 5
f <- 0.3
q <- 1
psi <- matrix(
c(
0.9, 0.05, 0.05,
0.05, 0.9, 0.05,
0.05, 0.05, 0.9
), nrow = 3, ncol = 3,
byrow = TRUE
)
par(mar = c(5,5,2,5))
plot(lambda[1, ], type = "l", col = "red", xlab = "Day", ylab = "Lambda (red)")
par(new = TRUE)
plot(psurv, type = "l", axes = F, xlab = NA, ylab = NA, col = "blue")
axis(side = 4)
mtext(side = 4, line = 3, 'Survival Probability (blue)')
## ---- echo=FALSE--------------------------------------------------------------
par(mar = pardefault)
## ---- fig.width=10, fig.height=8----------------------------------------------
M <- c(100, 100, 100)
Y <- c(0, 0, 0)
Z <- c(0, 0, 0)
mod <- make_MicroMoB(tmax = tmax, p = 3)
setup_mosquito_RM(mod, stochastic = FALSE, f = f, q = q, eip = EIP, p = psurv, psi = psi, M = M, Y = Y, Z = Z)
setup_aqua_trace(model = mod, lambda = lambda, stochastic = FALSE)
det_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), patch = 1:3, value = NaN)
det_out <- det_out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(det_out, day)
# run it
while(get_tnow(mod) <= tmax) {
mod$mosquito$kappa <- rep(0.05, 3)
step_aqua(model = mod)
step_mosquitoes(model = mod)
det_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
det_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
det_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
mod$global$tnow <- mod$global$tnow + 1L
}
ggplot(det_out) +
geom_line(aes(x = day, y = value, color = state)) +
facet_grid(patch ~ .)
## ---- fig.width=10, fig.height=8----------------------------------------------
sto_out <- mclapply(X = 1:10, FUN = function(runid) {
mod <- make_MicroMoB(tmax = tmax, p = 3)
setup_mosquito_RM(mod, stochastic = TRUE, f = f, q = q, eip = EIP, p = psurv, psi = psi, M = M, Y = Y, Z = Z)
setup_aqua_trace(model = mod, lambda = lambda, stochastic = TRUE)
out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), patch = 1:3, value = NaN)
out <- out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(out, day)
while(get_tnow(mod) <= tmax) {
mod$mosquito$kappa <- rep(0.05, 3)
step_aqua(model = mod)
step_mosquitoes(model = mod)
out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
mod$global$tnow <- mod$global$tnow + 1L
}
out[, 'run' := as.integer(runid)]
return(out)
})
sto_out <- data.table::rbindlist(sto_out)
ggplot(sto_out) +
geom_line(aes(x = day, y = value, color = state, group = interaction(run, state)), alpha = 0.35) +
facet_grid(patch ~ .)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.