ABC rejection algorithm

Below, you can find an example of how to code a ABC rejection algorithm. 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.

initState <- c(S = 250, E = 0, I = 4, T = 0, L = 30, Inc = 0)

my_ABCAlgorithm <- function(N, epsilon, sumStats, distanceAbc,
                            fitmodel, initState, data) {
  # set up empty matrix to store results
  results <- matrix(nrow = 0, ncol = 6)

  # initialise with i=0
  i <- 0

  # while the length of the accepted values (result) is less than the desired
  # length (N)
  while (i < N) {
    # - draw a new theta from prior distributions

    d_lat <- # INSERT HERE:
    d_inf <- # INSERT HERE:

    theta <- c(
      R_0 = 2, D_lat = d_lat, D_inf = d_inf, alpha = 0.9, D_imm = 13, rho = 0.85
    )

    # use computeDistanceAbc to calculate a distance between the model
    # and data
    dist <- computeDistanceAbc(
      sumStats = sumStats,
      distanceAbc = distanceAbc,
      fitmodel = fitmodel,
      theta = theta,
      initState = initState,
      data = data
    )

    ## if the model distance is within the epsilon window

    if (dist <= epsilon) {
      # store the accepted parameter values
      # INSERT HERE:
    }

    # update i (dimension of results store)
    i <- # INSERT HERE:
  }
  # return the accepted values
  return(results)
}

results <- my_ABCAlgorithm(
  N = 10, epsilon = 5,
  sumStats = list(ssMax, ssSize),
  distanceAbc = ssMeanRelDistance,
  fitmodel = seitlStoch,
  initState = initState,
  data = fluTdc1971
)

head(results)

hist(results[, 2])
hist(results[, 3])

If you run into any problems, have a look at our solution. Otherwise you can return to the ABC session.

Partial example: ABC rejection algorithm

ABC rejection algorithm