Bayesian Statistics

\CustomizeMathJax

8.4 Gibbs sampling

Recall that our parameter $θ$ may really be a vector of parameters, as in Section 4.5 where we considered the model $M_{(μ, σ)} \sim N (μ, σ^{2})$ with $θ = (μ, σ) \in R \times (0, \infty)$ . In this section we introduce a technique for handling models with many parameters $θ = (θ_{1}, \dots, θ_{d}) \in R^{d}$ . For those of you taking MAS364, this section is non-examinable and is included for interest only. For those of you taking MAS61006, this section is on syllabus and will feed into your work next semester.

For $d = 1$ or $d = 2$ the MH algorithm as described in Section 8.3 is effective. For much large $d$ , what tends to happen is that it takes a very long time for the sequence $(y_{m})$ generated by the MH algorithm to explore the parameter space $Π \subseteq R^{d}$ , which means that it takes longer to converge to (the distribution of) $Θ |_{{X = x}}$ . This happens simply because there is a lot of space to explore inside $R^{d}$ when $d$ is large; this phenomenon is often known as the curse of dimensionality.

Example 8.4.1 Imagine you have parked your car inside a multi-story car park and then forgotten where you’ve parked it. If you know which level your car is on then you will only have to search on that level, and you will find your car in minutes. If you don’t know which level, it will take you much longer. The first case is exploring $d = 2$ , the second is $d = 3$ . Actually, to make this example properly match the difference between $d = 2$ and $d = 3$ , the car park should have the same number of floors as there are parking spaces along the side-length of each floor! The problem only gets worse in $d \geq 3$ .

One strategy for working around this problem is to update the parameters $(θ_{1}, \dots, θ_{d})$ one at a time. That is, we would first change $θ_{1}$ while keeping $(θ_{2}, \dots, θ_{d})$ fixed, then we would change $θ_{2}$ while keeping $(θ_{1}, θ_{3}, \dots, θ_{d})$ fixed, and so on. After updating $θ_{d}$ we would then return to $θ_{1}$ . It is helpful to introduce some notation for this: we write $θ_{- i} = (θ_{1}, \dots, θ_{i - 1}, θ_{i + 1}, \dots, θ_{d})$ for the vector $θ$ with the $θ_{i}$ term removed. We use the same notation for the random vector $Θ$ e.g. $Θ_{- i}$ .

Remark 8.4.2 In reality, instead of updating the $θ_{i}$ one-by-one, it is common to update the parameters in small batches. For example we might update $(θ_{1}, \dots, θ_{4})$ in one step, then $(θ_{5}, \dots, θ_{8})$ in the next step, and so on. It is helpful to put related parameters, with values that might strongly influence each other, within the same batch.

When we update the parameters in turn, a common choice of proposal distribution is to set $Q \sim Θ_{i} |_{{Θ_{- i} = θ_{- i}, X = x}}$ in the update for $θ_{i}$ , where $θ_{- i}$ are the values obtained from the previous update. This choice of proposal has the effect that, from (8.1), we end up with $α = 1$ and all proposals are then accepted. When using proposals of this form, the MH algorithm is usually known as the Gibbs sampler.

Definition 8.4.3 The distributions of $Θ_{i} |_{{Θ_{- i} = θ_{- i}, X = x}}$ , for $i = 1, \dots, d$ , are known as the full conditional distributions. From Lemma 1.6.1 the $i^{t h}$ full conditional has p.d.f. given by

$\begin{aligned} f_{Θ_{i} |_{{Θ_{- i} = θ_{- i}, X = x}}} (θ_{i}) & = \frac{f_{Θ |_{{X = x}}} (θ)}{\int_{R^{d - 1}} f_{Θ |_{{X = x}}} (θ_{1}, \dots, θ_{i - 1}, θ_{i}, θ_{i + 1}, \dots, θ_{d}) d θ_{- i}} \\ (8.12) & \propto f_{Θ |_{{X = x}}} (θ) \end{aligned}$ We can calculate $f_{Θ |_{{X = x}}}$ from Theorems 2.4.1/3.1.2, which provides a strategy for calculating (8.12) analytically. Note that $\propto$ in (8.12) treats $θ_{- i}$ and $x$ as constants, and the only variable is $θ_{i}$ .

Remark 8.4.4 The notation $Θ_{i} |_{{Θ_{- i} = θ_{- i}, X = x}}$ for the full conditionals is a bit unwieldy. In Bayesian shorthand we would write simply $θ_{i} | θ_{- i}, x$ which is much neater.

The Gibbs sampler that results from these strategies is as follows.

1. Choose an initial point $y_{0} = (θ_{1}^{(0)}, \dots, θ_{d}^{(0)}) \in Π$ .
2. For each $i = 1, \dots, d$ , sample $\tilde{y}$ from $Θ_{i} |_{{Θ_{- i} = θ_{- i}^{(m)}, X = x}}$ and set

$y_{m + 1} = (θ_{1}^{(m)}, \dots, θ_{i - 1}^{(m)}, \tilde{y}, θ_{i + 1}^{(m)}, \dots, θ_{d}^{(m)}) .$

Note that we increment the value of $m$ each time that we increment $i$ . When reach $i = d$ , return to $i = 1$ and repeat.

Repeat this step until $m$ is large enough that values taken by the $y_{m}$ are no longer affected by the choice of $y_{0}$ .
3. The final value of $y_{m}$ is now a sample of $Θ |_{{X = x}}$ .

Note that we need to take samples from the full conditionals $Θ_{i} |_{{Θ_{- i} = θ_{- i}^{(m)}, X = x}}$ in step 2. This isn’t always possible, and the Gibbs sampler is only helpful if we can do that. It is often used in cases where the full conditionals turn out to be named distributions, or nearly one as in Example 8.4.5 below.

For the MH algorithm we had Theorem 8.2.2 to tell us that $y_{m}$ was (approximately) a sample of $Θ |_{{X = x}}$ , justified by the discussion in Section 8.2.3. It is possible to make similar arguments for the Gibbs algorithm above, but we won’t include them in our course.

If the full conditionals can’t be easily sampled from, then one strategy is to use the MH algorithm (run inside of step 2 above) to obtain samples of $Θ_{- i} |_{{X = x}}$ . This technique is known as Metropolis-within-Gibbs. In practice, once the parameters are divided up into batches, as described in Remark 8.4.2, some batches may be amenable to Gibbs sampling, whilst others may require Metropolis-within-Gibbs. The details will depend on the model. Trial and error is often required to find the best combination of techniques. We won’t try to write down algorithms of that complexity within these notes, but you should hopefully end the course with an understanding of how (and why) each piece of an algorithm like that would work.

Example 8.4.5 This example comes from Sections 1.1.1/7.5.3/8.6.2 of the book ‘Bayesian Approach to Intrepreting Archaeological Data’ by Buck et al (1996). The data comes from radiocarbon dating, and is a vector $(x_{1}, x_{2}, \dots, x_{n})$ of estimated ages obtained (via carbon dating) from $n$ different objects. We write $θ_{i}$ for the true age of object $i$ , which is unknown. Our model for the age of each object is $x_{i} \sim N (θ_{i}, v_{i})$ and we assume that the estimation errors are independent, for each $i$ . For simplicity we will assume that the $v_{i}$ are known parameters, so our model family has $n$ parameters $θ = (θ_{1}, \dots, θ_{n})$ . We thus have the model family

$M_{θ} = N (θ_{1}, v_{1}) \otimes \dots \otimes N (θ_{n}, v_{n}) .$

From the historical context of the objects, it is known that $θ_{1} < θ_{2} < \dots < θ_{n}$ , so we condition our model $M_{θ}$ on this event. We can use Exercise 1.8 to do this conditioning, resulting in a new model family $M_{θ}^{'}$ given by

$\begin{aligned} f_{M_{θ}^{'}} (x) & = {\begin{cases} \frac{f_{M_{θ}} (x)}{P [N (θ_{1}, v_{1}) < N (θ_{2}, v_{2}) < \dots < N (θ_{n}, v_{n})]} & for θ_{1} < θ_{2} < \dots < θ_{n} \\ 0 & otherwise \end{cases} \\ \propto {\begin{cases} \prod_{i = 1}^{n} f_{N (θ_{i}, v_{i})} (x_{i}) & for θ_{1} < θ_{2} < \dots < θ_{n} \\ 0 & otherwise \end{cases} \\ \propto {\begin{cases} \exp (- \frac{1}{2} \sum_{i = 1}^{n} \frac{(θ_{i} - x_{i})^{2}}{v_{i}}) & for θ_{1} < θ_{2} < \dots < θ_{n} \\ 0 & otherwise. \end{cases} \end{aligned}$ We use the Bayesian model $(X, Θ)$ with model family $M_{θ}^{'}$ and the improper prior

$f_{Θ} (θ) = {\begin{cases} 1 & for 0 < θ_{1} < θ_{2} < \dots < θ_{n}, \\ 0 & otherwise. \end{cases}$

By Theorem 3.1.2 we obtain that the posterior distribution has p.d.f.

$\begin{matrix} (8.13) & f_{Θ |_{{X = x}}} (θ) \propto {\begin{cases} \exp (- \frac{1}{2} \sum_{i = 1}^{n} \frac{(θ_{i} - x_{i})^{2}}{v_{i}}) & for 0 < θ_{1} < θ_{2} < \dots < θ_{n}, \\ 0 & otherwise. \end{cases} \end{matrix}$

This is the same density as $M_{θ}^{'}$ , except now we treat $θ$ rather than $x$ as the variable. The density is symmetric in $x$ and $θ$ so we already know this distribution. It is the distribution of $θ \sim N (x_{1}, v_{1}) \otimes \dots \otimes N (x_{n}, v_{n})$ conditioned on the event $0 < θ_{1} < θ_{2} < \dots < θ_{n}$ . One way to simulate samples of this distribution is via rejection sampling: simulate $θ \sim N (x_{1}, v_{1}) \otimes \dots \otimes N (x_{n}, v_{n})$ and reject the sample $θ$ until it satisfies $0 < θ_{1} < θ_{2} < \dots < θ_{n}$ . The trouble is that unless $n$ is small, we will mostly end up rejecting the samples because the condition we have imposed is an unlikely one.

From (8.13) and (8.12) we have full conditionals given by

$\begin{aligned} f_{Θ_{i} |_{{Θ_{- i} = θ_{- i}, X = x}}} (θ_{i}) & \propto {\begin{cases} \exp (- \frac{1}{2} \sum_{j = 1}^{n} \frac{(θ_{j} - x_{i})^{2}}{v_{i}}) & for θ_{i} \in (θ_{i - 1}, θ_{i + 1}), \\ 0 & otherwise \end{cases} \\ \propto {\begin{cases} \exp (- \frac{1}{2} \frac{(θ_{i} - x_{i})^{2}}{v_{i}}) & for θ_{i} \in (θ_{i - 1}, θ_{i + 1}), \\ 0 & otherwise \end{cases} \end{aligned}$ (where we set $θ_{0} = 0$ and $θ_{n + 1} = \infty$ to make convenient notation). Note that $θ_{i}$ is the only variable here, and the second line follows because $θ_{- i}$ and $x$ are treated as constants. We recognize this full conditional distribution as that of $θ_{i} \sim N (x_{i}, v_{i})$ conditioned on the event $θ_{i} \in (θ_{i - 1}, θ_{i + 1})$ . These full conditionals are much easier to sample from: we use rejection sampling, sample $θ_{i} \sim N (x_{i}, v_{i})$ and reject until we obtain a sample for which $θ_{i} \in (θ_{i - 1}, θ_{i + 1})$ . Hence, in this situation we have all the necessary ingredients to use a Gibbs sampler.