In Lecture 2 we discussed about the modelling process in general and what steps does one have to follow to perform Bayesian inference. In this lab session we look at how to implement an MCMC algorithm to infer the parameters that govern the incubation period distribution that was discussed in the first part of the lecture.
Suppose that we have data on incubation periods \(y = (y_1, \ldots, y_n)\). We assume that the data are independent draws from a Gamma distribution with shape \(\alpha\) and rate \(\beta\), i.e. \[y_i \sim Ga(\alpha, \beta), \qquad \mbox{(i.i.d.)} \quad i=1, \ldots, n.\] where the Gamma distribution has probability density function: \[ f(y|\alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{\alpha-1} \exp\{-\beta y\}, \qquad \alpha > 0, \quad \beta > 0, \qquad y > 0. \]
The goal is to make (sampling-based) Bayesian inference for the parameters \(\alpha\) and \(\beta\).
Assume that \(\alpha\) and \(\beta\) are a priori independent and assign Gamma distributions with parameters \(\lambda_{\alpha}\) and \(\nu_\alpha\), and \(\lambda_{\beta}\) and \(\nu_\beta\), respectively. Write down the posterior density of interest, i.e. \[ \pi(\alpha, \beta | y) \propto \pi(y|\alpha, \beta) \pi(\alpha) \pi(\beta) \]
where \(\pi(y|\alpha, \beta)\) denotes the likelihood function.
Derive the densities of the posterior distribution of the parameters \(\alpha\) and \(\beta\) up to proportionality, i.e. \(\pi(\alpha|\beta, y)\) and \(\pi(\beta|\alpha, y)\).
Is any of these densities of a known form (i.e. the density of a standard distribution)?
Write an MCMC algorithm (from scratch!) which draws samples from the joint posterior distribution \(\pi(\alpha, \beta|y)\).
A pseudo-code of this MCMC algorithm is given below:
- Choose initial values for \(\alpha\) and \(\beta\);
- Update \(\alpha\) by sampling from \(\pi(\alpha|\beta, y)\) using a Gaussian random-walk;
- Update \(\beta\) by sampling from \(\pi(\beta|\alpha, y)\) directly.
- Goto to Step (b)
Look at the trace plots and convince yourself that the chain has reached stationarity and is mixing well.
Test your algorithm is producing sensible results.
One way to do this is to first simulate a dataset from a Gamma
distribution with some specific values for \(\alpha\) and \(\beta\), eg.
y <- rgamma(100, 4, 2) wil simulate a vector of 100
draws from a Ga\((4, 2)\).
Then use your MCMC algorithm (in Exercise 3) to sample from the joint posterior distribution of \(\alpha\) and \(\beta\) assuming vague priors, e.g. \(\lambda_\alpha = \lambda_\beta = 1\) and \(\nu_\alpha = \nu_\beta = 10^{-3}\).
In principle, your marginal posterior distributions (\(\pi(\alpha|y)\) and \(\pi(\beta|y)\)) should be centered around the values you have chosen to simulate the data from (e.g. \(4\) and \(2\) in this case).
If you have convinced yourself that your MCMC algorithm is doing what is supposed to be doing, then use it to fit the model to the Campylobacter data (24 observations in total) from Evans et al. 1996) by
y <- c(rep(2,2), rep(3, 6), rep(4, 11), rep(5, 3), rep(7,7))Does the model offer a good fit? How would you go assessing that?
Now that you have got a working algorithm, can you improve its mixing?
One way to do this is to update both \(\alpha\) and \(\beta\) at the same time (in a block). This can be done, for example, and as discussed in the lecture, by proposing values drawn from some (bivariate distribution) with density \(q(\alpha, \beta)\) and accept/reject according to the Metropolis-Hastings ratio. There are many choices here, but not all of them will lead to an efficient sampler.
Here are a few choices to try:
optim to find the values of \(\alpha\) and \(\beta\) that maximise the posterior
density).