Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup--------------------------------------------------------------------
library(MicroMoB)
library(ggplot2)
library(data.table)
library(parallel)
## -----------------------------------------------------------------------------
h <- 0.025
r <- 1/200
b <- 0.55
EIR <- -log(1 - h) / b
n <- 1
tmax <- 1e3
MOI_init <- matrix(data = c(1e5, rep(0, 1e2)), nrow = 101, ncol = n)
mod <- make_MicroMoB(tmax = tmax, p = 1)
setup_humans_MOI(model = mod, stochastic = TRUE, theta = matrix(1, nrow = n, ncol = 1), H = colSums(MOI_init), MOI = MOI_init, r = r, b = b)
human_out <- data.table::CJ(day = 1:tmax, MOI = 0:(nrow(MOI_init)-1), value = NaN)
data.table::setkey(human_out, day)
while (get_tnow(mod) <= mod$global$tmax) {
mod$human$EIR <- EIR
step_humans(model = mod)
human_out[day==get_tnow(mod), value := as.vector(mod$human$MOI)]
mod$global$tnow <- mod$global$tnow + 1L
}
## -----------------------------------------------------------------------------
weighted.mean(x = 0:100, w = human_out[day == tmax, value])
## -----------------------------------------------------------------------------
ggplot(data = human_out) +
geom_line(aes(x = day, y = value, group = MOI, color = MOI))
## -----------------------------------------------------------------------------
human_final <- human_out[day == tmax, ]
human_final[, 'theoretical' := dpois(x = MOI, lambda = h/r)]
human_final[, 'empirical' := value / sum(value)]
ggplot(human_final, aes(MOI, empirical)) +
geom_bar(stat = 'identity') +
geom_line(aes(x = MOI, y = theoretical), color = "blue")
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.