Background

The purpose of this lab session is to illustrate the general idea of Approximate Bayesian Computation algorithms with application to the epidemic models that we have been looking at in the lectures and the lab sessions.

In particular, we will look at the data from the Lab session 4 on Asian influenza taken from Longini & Koopman (1982). Recall that there is a community which consists of 42 households each of them of size 3. ```{r}

cases <- c(0,1,2,3)
households <- c(29,9,2,2)
mat <- matrix(NA, nrow=4, 2)
mat[,1] <- cases
mat[,2] <- households
mat <- as.data.frame(mat)
colnames(mat) <- c("No.of.Cases", "No.of.Households")
mat
##   No.of.Cases No.of.Households
## 1           0               29
## 2           1                9
## 3           2                2
## 4           3                2

We want to fit the model that we fitted before to the same data in Lab session 4 using ABC. Although we can compute the likelihood of the observed data given the parameters, we will pretend that we don’t for the purpose of illustration.

Simulation

The function below simulates data from the model

sim.final.size.data <- function(alpha, p) {
  
  # compute final size probabilities
  p1 <- compute.marginal.prob(0,3,1,alpha, p)
  p2 <- compute.marginal.prob(1,3,1,alpha, p)
  p3 <- compute.marginal.prob(2,3,1,alpha, p)
  p4 <- compute.marginal.prob(3,3,1,alpha, p)

  # put them in a vector
  final.size.prob.vector <- c(p1, p2, p3, p4)
  
  # draw from a multinomial distribution
  sim.data <- rmultinom(1, size = 42, final.size.prob.vector)
  return(sim.data)
  
}

Exercises

Exercise 1

Assume that a priori \[p \sim \mbox{Beta}(1, 1) \qquad \mbox{ and } \qquad \alpha \sim \mbox{Gamma}(1, 1).\]

Write an algorithm to draw exact samples from the posterior distribution \[\pi(\alpha, p| y),\] where \(y\) denotes the observed data, using Algorithm ER2 (ie exact rejection sampling). You may find useful to look at the pseudo-code in the notes from Lecture 5.2

Exercise 2

Assuming the same priors as above, write an algorithm to draw samples from an approximation to the posterior distribution. In other words, write some code to implemnt an ABC algorithm.

Have a think as to what might be appropriate to use for comparing how “close” the simulated and the observed data are.

Exercise 3

Compare how good is your approximate posterior distribution to the true one (i.e. the one you draw samples from using MCMC in Lab Session 4).

Exercise 4

Use your ER2 and ABC algorithm in Exercises (2 and 3) assuming a vague prior for \(\alpha\), e.g. \(\alpha \sim Ga(1,0.001)\) and comment on the efficiency of both algorithms.