Lab 7/Homework 10

Due Wednesday, December 3, 11:59pm

A common model in epidemiology is a Susceptible-Infected-Recovered (SIR) model. These are usually written as deterministic, but a simple modification makes it into a Markov chain, which is slightly more realistic.

In this assignment, you will simulate from an SIR Markov chain.

As we discussed in class, to define a Markov chain we need to define the state space and the transition probabilities. For our SIR model, the state space is going to be the set of triples \((S, I, R)\), such that \(S, I, R \in \{0, \ldots, N\}\), and \(S + I + R = N\). The idea is that we have a total population of \(N\) individuals, with \(S\) of those being susceptible to the disease, \(I\) of those being currently infected, and \(R\) being recovered and no longer susceptible to the disease.

For the transition probabilities, from one time step to the next we can either move one individual from the susceptible compartment to the infected compartment (\((S, I, R) \rightarrow (S - 1, I + 1, R)\)) or we can move one individual from the infected compartment to the recovered compartment (\((S, I, R) \rightarrow (S, I - 1, R + 1)\)). If \(S + I > 0\), the probabilities are \[\begin{align} P[(S, I, R) \rightarrow (S - 1, I + 1, R)] &= \beta \frac{I}{N} S / C \label{eq:pi}\\ P[(S, I, R) \rightarrow (S, I - 1, R + 1)] &= \gamma I / C \label{eq:pr}\\ C &= \beta \frac{I}{N} S + \gamma I, \end{align}\] where \(\beta\) and \(\gamma\) are parameters of the model. If \(S + I = 0\), the transition probability is: \[\begin{align} P[(S, I, R) \rightarrow (S, I, R)] = 1.\label{eq:pab} \end{align}\] All other transitions have probability 0.

The idea is that \(\gamma\) controls the recovery rate (large values of \(\gamma\) mean that individuals leave the “infected” compartment quickly), and \(\beta\) controls the infection rate. The state \((0,0,N)\) is an “absorbing” state: once everyone has recovered, the chain stays at \((0,0,N)\) forever.

What you should do:

Write a function that takes a current state (that is, a triple \((S, I, R)\)), values for the parameters \(\beta\) and \(\gamma\), and simulates one step in the Markov chain.

## current_state should be a list with elements S, I, and R.  parameters
## should be a list containing beta and gamma, the model parameters
take_one_step <- function(current_state, parameters) {

}

That is, your function will compute the probabilities in (\(\ref{eq:pi}\)-\(\ref{eq:pab}\)), and will return either \((S - 1, I + 1, R)\) or \((S, I - 1, R + 1)\) with the probabilities you computed.

Write a function that will simulate some fixed number of step (say, 1000), from the SIR Markov chain. The function should take a starting state (for example, \((99, 1, 0)\)) and should record the state of the chain across all the steps.

## starting_state should be a list with elements S, I, and R.
## parameters should be a list containing beta and gamma, the model
## parameters. n_steps should be the number of steps to take. It will
## probably be easiest to deal with if the output is a list of
## vectors, one for S, one for I, and one for R.
simulate_sir <- function(starting_state, parameters n_steps) {
    
}

Run your chain and either plot or print out some summary of how the chain evolved.
Your simulation takes fixed values of \(\beta\) and \(\gamma\). Suppose you had data from a real epidemic and wanted to model the real epidemic using this SIR model.

draw_params_from_prior <- function() {
    beta <- runif(0, 10)
    gamma <- runif(0, 10)
    return(list(beta = beta, gamma = gamma))
}
## sir_path should be in whatever format the output of simulate_sir
## is. The output should be a vector of summary statistics.
summary_statistic <- function(sir_path) {

}

Submission parameters

Submit two files:

A pdf or html containing your plots and answers to the questions.
An Rmd or R file containing the code you used.