## 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)