library(deSolve)
library(magrittr)
library(reshape2)
library(functional)
library(ggplot2)
## StepSir: Given Y_{t} and parameters of transition, return delta(Y_{t+1})
## times argument is not used, but necessary?
## Output a list corresponding to delta(Y_{t+1})
StepSir <- function(times, Yt, parms) {
with(as.list(c(Yt, parms)), {
## delta of each compartment
dSusceptible <- birth - beta*Infected*Susceptible - death*Susceptible
dInfected <- beta*Infected*Susceptible - recovery*Infected - death*Infected
dRecovered <- recovery*Infected - death*Recovered
## delta(Y_{t + 1})
c(dSusceptible, dInfected, dRecovered)
}) %>% list
}
## Partial application to obtain a function that takes parameters only
SimulateSir <- Curry(FUN = ode,
func = StepSir,
y = c(Susceptible = 0.9, Infected = 0.1, Recovered = 0),
times = seq(from = 0, to = 1000, by = 1))
## Function to run simulation and plot
RunSimulation <- function(parms) {
## Run simulation
dfMelt <- SimulateSir(parms = parms) %>%
as.data.frame %>%
## Melt
melt(data = . ,
id.vars = "time",
variable.name = "variable",
value.name = "value")
## Return ggplot2 plot object
ggplot(data = dfMelt, mapping = aes(x = time, y = value, color = variable)) +
layer(geom = "line") +
theme_bw() + theme(legend.key = element_blank())
}
## Run simulations
## Balanced birth and death
RunSimulation(parms = c(beta = 0.1, recovery = 0.005, death = 0.001, birth = 0.001)) %>% plot
## No birth or death
RunSimulation(parms = c(beta = 0.1, recovery = 0.005, death = 0.000, birth = 0.000)) %>% plot
## Faster recovery
RunSimulation(parms = c(beta = 0.1, recovery = 0.01, death = 0.000, birth = 0.000)) %>% plot
