Skip to content

Exercise #3: Markov model

We will now look at implementing a Markov model, specifically a simple SIR (susceptible-infectious-recovered) compartmental model for an influenza epidemic, using a classic data set of an outbreak of Russian influenza at a boarding school.

Model

The model is described in three parts: a parameter model, an initial model, and a transition model. Time is indexed by , in days. The state consists of variables , , and , giving counts of the number of susceptible, infectious, and recovered individuals, respectively.

Parameter model

The parameter model is: where is a rate of interaction in the population, the probability of infection when a susceptible individual interacts with an infectious individual, and the daily recovery probability.

Initial model

The initial model for time depends on the data set. For the data set introduced below, it is as follows:

Transition model

The transition model for time is: where is the number of interactions between infectious and susceptible individuals, the number of newly infected individuals, and the number of newly recovered individuals. Population counts are then updated:

Implementation

To specify this model in Birch, we again need to create a class that inherits from Model. This time, however, the standard library provides the more-specific class MarkovModel that will do some of the work for us. MarkovModel itself inherits from Model, we just need to provide it with specific implementations of the parameter, initial and transition models.

As well as easing implementation, the other advantage of MarkovModel is that it reveals something about the structure of the model, which may be useful to enable, or optimize, specific inference methods. It is expected that more classes for specific model structures, and more methods to make use of them, will be available in future.

Exercise

Create a file bi/SIRModel.bi, add it to META.json, and enter the following:

/**
 * SIR (susceptible-infectious-recovered) model.
 */
class SIRModel = MarkovModel<SIRState,SIRParameter>;

MarkovModel is a generic class. It takes the name of two other classes between the angle brackets—in this case SIRState and SIRParameter—which it uses internally. We will create these classes below. Otherwise, this code is just establishing the name SIRModel as an alias for MarkovModel<SIRState,SIRParameter>: two names for the same type.

We will start with the SIRParameter class. This is the parameter model. It must inherit from Parameter. The Parameter class is very similar to the Model class, but it provides a parameter member fiber that must be overriden to define the parameter model.

Exercise

Create a file bi/SIRParameter.bi, add it to META.json, and enter the following:

/**
 * SIR model parameters.
 */
class SIRParameter < Parameter {
  /**
   * Interaction rate.
   */
  λ:Random<Real>;

  /**
   * Infection probability.
   */
  δ:Random<Real>;

  /**
   * Recovery probability.
   */
  γ:Random<Real>;

  fiber parameter() -> Real {
    λ <- 10.0;
    δ ~ Beta(2.0, 2.0);
    γ ~ Beta(2.0, 2.0);
  }

  function input(reader:Reader) {
    λ <- reader.getReal("λ");
    δ <- reader.getReal("δ");
    γ <- reader.getReal("γ");
  }

  function output(writer:Writer) {
    writer.setReal("λ", λ);
    writer.setReal("δ", δ);
    writer.setReal("γ", γ);
  }
}

Recall the typical structure of this class from the Bayesian linear regression example: random variables as member variables, the simulate member fiber, the input and output member functions.

The SIRState class must inherit from State. The State class is very similar to the Model class, but it provides two separate member fibers initial and transition that must be overridden to define the initial and transition models, respectively.

Exercise

Create a file bi/SIRState.bi, add it to META.json, and enter the following:

/**
 * SIR model state.
 */
class SIRState < State {
  /**
   * Number of susceptible-infectious interactions.
   */
  τ:Random<Integer>;

  /**
   * Newly infected population.
   */
  Δi:Random<Integer>;

  /**
   * Newly recovered population.
   */
  Δr:Random<Integer>;

  /**
   * Susceptible population.
   */
  s:Random<Integer>;

  /**
   * Infectious population.
   */
  i:Random<Integer>;

  /**
   * Recovered population.
   */
  r:Random<Integer>;

  fiber initial(θ:SIRParameter) -> Real {
    //
  }

  fiber transition(x:SIRState, θ:SIRParameter) -> Real {
    τ ~ Binomial(x.s, 1.0 - exp(-θ.λ*x.i/(x.s + x.i + x.r)));
    Δi ~ Binomial(τ, θ.δ);
    Δr ~ Binomial(x.i, θ.γ);

    s ~ Delta(x.s - Δi);
    i ~ Delta(x.i + Δi - Δr);
    r ~ Delta(x.r + Δr);
  }

  function input(reader:Reader) {
    Δi <- reader.getInteger("Δi");
    Δr <- reader.getInteger("Δr");
    s <- reader.getInteger("s");
    i <- reader.getInteger("i");
    r <- reader.getInteger("r");
  }

  function output(writer:Writer) {
    writer.setInteger("Δi", Δi);
    writer.setInteger("Δr", Δr);
    writer.setInteger("s", s);
    writer.setInteger("i", i);
    writer.setInteger("r", r);
  }
}

Notice the two member fibers in the above code:

fiber initial(θ:SIRParameter) -> Real;
fiber transition(x:SIRState, θ:SIRParameter) -> Real;

As suggested by their names and parameters:

  • the first is for the initial model, providing the parameters as the θ argument,
  • the second is for the transition model, providing both the previous state as the x argument and the parameters as the θ argument.

The initial model is empty by choice. While we could implement the initial model described above, it is specific to the data set that we will use. We would prefer not to hardcode it, in order that we might reuse this model for other data sets. Instead, we have elected to include the initial state in the input file that we will set up below.

The transition model uses Delta distributions rather than simple assignment statements. This is because we have declared the variables s, i and r to be of type Random<Integer>, not just Integer. Using the Delta distributions ensures proper handling of the two possible situations for each variable: that it already has a value—in which case we observe that value—or that its value is missing—in which case we simulate it. Indeed, the code, as written, handles all cases, according to what is provided in the input file.

Exercise

Build the project with

birch build

Data

We will use a data set of the outbreak of Russian influenza at a boy's boarding school in northern England1.

Exercise

Download the data set here and place it in your project's input/ directory as input/russian_influenza.json.

Add the file to META.json under manifest.data.

Have a look at the contents of the file in a text editor. It contains an array of states. The first state sets the values of all state variables, while for subsequent states it sets only , the number of infectious individuals.

Inference

We can now run the model.

Exercise

Sample from the posterior distribution with

birch sample \
  --model SIRModel \
  --input-file input/russian_influenza.json \
  --output-file output/russian_influenza.json \
  --ncheckpoints 14 \
  --nparticles 100 \
  --nsamples 10

The new command-line option --ncheckpoints gives the number of checkpoints for which to run. In the case of a model that inherits from MarkovModel, as here, this is interpreted as the number of states. The numbers that appear in the terminal are simply ticking off these checkpoints as the sampler proceeds.

Unlike the previous example, there is no exact analytical solution for this model that can be computed in reasonable time. A particle filter is used for inference instead. The new command-line option --nparticles gives the number of particles to use in the particle filter.

The delayed sampling2 mechanism in Birch automatically applies some analytical optimizations in this case too. These include marginalizing out the parameters and , and enumerating sums and differences of binomials.

As before, you can inspect the results of the inference in output/russian_influenza.json, and perhaps plot them in a package such as MATLAB, R, or Julia. Be aware that the output here is an importance sample. Each sample is assigned a weight, the logarithm of which is given by the associated weight element in the output file.


  1. Anonymous (1978). Influenza in a boarding school. British Medical Journal. 1:587. 

  2. L.M. Murray, D. Lundén, J. Kudlicka, D. Broman and T.B. Schön (2018). Delayed Sampling and Automatic Rao–Blackwellization of Probabilistic Programs. In Proceedings of the 21st International Conference on Artificial Intelligence and Statistics (AISTATS) 2018, Lanzarote, Spain.