## Headache example: discrete model of a non-infectious process ## Last modified: Jan 24 2018 ## Let's set some initial values and parameters N_t <- 500 # Number of people without a headache at time zero P_t <- 0 # Number of people with a headache at time zero incidence <- 0.05 # Probability of going from {No Headache} to {Headache} recovery <- 0.9 # Probability of going from {Headache} to {No Headache} ## Again, we make a data container dat <- NULL ## Make a sequence of time steps timesteps <- 1:100 ## Or, equivalently, seq(1, 100, 1) ## Loop through our sequence of time steps for (t in timesteps){ ## Calculate the number of unaffected at time t + 1 N_t1 <- N_t - round(incidence*N_t) + round(recovery*P_t) ## Calculate the number of affected at time t + 1 P_t1 <- P_t + round(incidence*N_t) - round(recovery*P_t) ## Bind these numbers as a new row (rbind) in our data container dat <- rbind(dat, c(N_t1, P_t1)) ## Update the unaffected and affect for the next time step N_t <- N_t1 P_t <- P_t1 } ## Plot it matplot(dat, xlab="time", ylab="number of people", type="l", lty = 1) legend("topright", col = 1:2, legend = c("no headache","headache"), lwd=1)