Classes will be held 13:30-14:45 on Mondays and Wednesdays, in
Harrelson Hall, Room 368.
MATLAB Code for various parts of the course
Deterministic simulation code:
sir_det.m code that simulates SIR model
in a closed population and produces time course plots and phase plot.
sir_det_phase.m produces phase plot
of SIR model in closed population, using a number of different initial
conditions.
sir_demography_det.m code that
simulates the SIR model with demography. (Artifically high birth rate is
chosen for clarity. See code for more discussion.)
Stochastic simulation code:
sir_closed.m code that does a single
realization of the SIR model in a closed population. (The only processes
considered are infection and recovery.)
sir_stoch.m code to
generate outbreak size histograms for the closed SIR model.
sir_stoch_next_event.m code that
produces outbreak size histograms using the next event simulation method.
Stochastic simulation is slow in MATLAB for all but the smallest
population sizes.
Here is some FORTRAN code that does the stochastic simulations much more
quickly:
sir_stoch.f. SIR model in closed
population. (It seems to run several orders of magnitude faster than
MATLAB.)
sir_stoch_demog.c SIR model
with demography. Carries out one realization, starting at the endemic
equilibrium. Code is designed to print out (S,I) at regular intervals.
sir_stoch_variability.c
Calculates variability in the SIR model with demography. Moments are
calculated numerically by producing a collection of 1000 realizations.
Branching process code:
branching_process_poisson.m.
Code that simulates the branching process with Poisson distributed numbers
of secondary infections (corresponds to each individual having a fixed
duration of infection, during which time they give rise to secondary
infections at a constant rate).
invasion_prob_poisson.m.
Code that repeatedly carries out the previous simulation, and from these
realizations estimates the probability of successful invasion of the
infection.
invasion_prob_poisson_graph.m. Code that estimates the invasion
probability for a number of values of R_0, and plots a graph.
invasion_prob_geometric.m.
Code that estimates invasion probability assuming that secondary
infections follow a geometric distribution. (This arises when infectious
periods are exponentially distributed, and infectious individuals give
rise to secondary infections at a constant rate over their infectious
period.) In this case, we have the simple formula for invasion
probability: its predicted value is outputted.
invasion_prob_geometric_graph.m. Code that plots invasion probability
as a function of R_0. Makes use of the fact that the sum of independent,
identically distributed, geometric random variables is distributed
according to the negative binomial distribution.
Chain Binomial Code:
max_likelihood_greenwood.m.
Calculates maximum likelihood estimate for the probability q in the
Greenwood model for the Providence measles data, based on outbreak sizes.
max_likelihood_reed_frost.m.
Calculates maximum likelihood estimate for the probability q in the
Reed-Frost model for the Providence measles data, based on outbreak sizes.
SIR model with fixed infectious period:
sir_delay.m. Compares the behavior of the
SIR model with fixed infectious period to that of the standard SIR
model (exponential infectious period). The latter is shown in red. Notice
how the epidemic takes off more quickly in the former case.
You might also see that the number of infectives can become negative
in the fixed period simulation. This is because MATLAB interpolates to
calculate S(t-tau) and I(t-tau). Numerical errors in that interpolation,
and in the discretization of the ODE, can lead to more people recovering
at time t than were infected at time t-tau.