ABC rejection algorithm

Here is a possible solution for the ABC rejection algorithm.

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 <- rgamma(1, shape = 16, rate = 8)
    d_inf <- rgamma(1, shape = 16, rate = 8)

    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
      results <- rbind(results, theta)
    }

    # update i (dimension of results store)
    i <- dim(results)[1]
  }
  # return the accepted values
  return(results)
}

You can copy and paste the function into your R session, and proceed from there.

Return to the ABC session.

Full example: ABC rejection algorithm

ABC rejection algorithm