Partial example: The particle filter

Below, you can find an example of how to code a particle filter. Some bits are left out for you to fill in (marked “INSERT HERE”). Each “INSERT HERE” statement requires one line of code. If you struggle, you can find a link to a solution below the function.

# This is a function that takes four parameters:
# - fitmodel: a fitmodel object
# - theta: named numeric vector. Values of the parameters for which the marginal log-likelihood is desired.
# - initState: named numeric vector. Initial values of the state variables.
# - data: data frame. Observation times and observed data.
# The function returns the value of the marginal log-likelihood
my_particleFilter <- function(fitmodel, theta, initState, data, nParticles) {
  ## Initialisation of the algorithm

  # Marginal log-likelihood is set to 0 and will be updated during the filtering steps
  margLogLike <- 0

  # Particle states can be stored in a list
  stateParticles <- # INSERT HERE

  # Weight: initially equal for all the particles
  # particle weight can be stored in a vector
  weightParticles <- # INSERT HERE

  # Initialise time variable
  currentTime <- 0

  ## Loop over observation times: resample, propagate, weight
  for (i in seq_len(nrow(data))) {
    # Extract next data point (must be a vector)
    dataPoint <- # INSERT HERE
    nextTime <- dataPoint["time"]

    # Resample particles according to their weights.
    # You can use the `sample` function of R (normalization of the weights is done in the function)
    indexResampled <- # INSERT HERE
    stateParticles <- stateParticles[indexResampled]

    ## Loop over particles: propagate and weight
    for (p in 1:nParticles) {
      # Extract current state of the particle
      currentStateParticle <- stateParticles[[p]]

      # Propagate the particle from current observation time
      # to the next one using the function `fitmodel$simulate`
      traj <- # INSERT HERE

      # Extract state of the model at next observation time.
      # Make sure that modelPoint is a vector
      modelPoint <- # INSERT HERE

      # Weight the particle with the likelihood of the observed
      # data point using the function `fitmodel$dPointObs`
      weightParticles[p] <- # INSERT HERE

      # Update state of the p particle
      stateParticles[[p]] <- modelPoint
    }

    # Increment time
    currentTime <- nextTime

    ## Increment the marginal log-likelihood
    # Add the log of the mean of the particles weights
    margLogLike <- margLogLike + # INSERT HERE
  }

  ## Return marginal log-likelihood
  return(margLogLike)
}

If you run into any problems, have a look at our solution, otherwise return to the pMCMC practical.

 

---

This web site and the material contained in it were originally created in support of an annual short course on Model Fitting and Inference for Infectious Disease Dynamics at the London School of Hygiene & Tropical Medicine. All material is under a MIT license. Please report any issues or suggestions for improvement on the corresponding GitHub issue tracker. We are always keen to hear about any uses of the material here, so please do get in touch using the Discussion board if you have any questions or ideas, or if you find the material here useful or use it in your own teaching.