Bayesian Statistics
Chapter 8 Computational methods
We noted at several points that conjugate pairs, from Chapter 4, do not provide enough flexibility for many practical situations. Instead, Bayesian statistics is heavily reliant on a family of computational techniques, introduced within this chapter. They generate samples from the posterior distribution and have a moderate computational cost. With simple model families of the kind we have worked with throughout the course it is reasonable to use desktop machines. More complex model families can require larger machines.
Throughout this chapter we assume the same setup as in Chapter 7, which we repeat here for convenience. We work with a discrete or absolutely continuous Bayesian model \((X,\Theta )\), where we have data \(x\) and posterior \(\Theta |_{\{X=x\}}\). We keep all of our usual notation: the parameter space is \(\Pi \), the model family is \((M_\theta )_{\theta \in \Pi }\), and the range of the model is \(R\). Note that \(M_\theta \) could have the form \(M_\theta \sim (Y_\theta )^{\otimes n}\) for some random variable \(Y_\theta \) with parameter \(\theta \), corresponding to \(n\) i.i.d. data points.
8.1 Approximate Bayesian computation \(\offsyl \)
In this section we describe a numerical method for calculating the posterior \(\Theta |_{\{X=x\}}\) that is based on rejection sampling. Recall that we used rejection sampling in Section 1.4 to prove Lemma 1.4.1, and also to give some intuition for our first examples of conditioning.
The algorithm we study here is known as Approximate Bayesian Computation, or ABC for short. We will describe it first for discrete data, in the situation where the prior (and, consequently, the posterior) are also discrete distributions. We haven’t studied this case in any of our previous chapters, so let us first introduce it here.
-
Definition 8.1.1 (Bayesian model with discrete parameters and discrete data) Take a prior with p.m.f. \(p_\Theta (\theta )\) and a discrete model family \((M_\theta )_{\theta \in \Pi }\) where \(\Pi \) is a finite or countable set. The Bayesian model \((X,\Theta )\) has the law
\[\P [X=x,\Theta =\theta ]=\P [M_\theta =x]p_\Theta (\theta ).\]
It is straightforward to sum over \(x\) and obtain the prior distribution \(\P [\Theta =\theta ]=p_\Theta (\theta )\), and also to sum over \(\theta \) and obtain the sampling distribution \(\P [X=x]=\sum _{\theta \in \Pi }\P [M_\theta =x]p_\Theta (\theta )\). For \(x\in R_X\) we have \(\P [X=x]>0\) and thus the posterior \(\Theta |_{\{X=x\}}\) is defined via Lemma 1.4.1. Also using Lemma 1.4.1, the conditional distribution \(X|_{\{\Theta =\theta \}}\) satisfies
\[\P [X|_{\{\Theta =\theta \}}=x] =\frac {\P [X=x,\Theta =\theta ]}{\P [\Theta =\theta ]} =\frac {\P [M_\theta =x]p_\Theta (\theta )}{p_{\Theta }(\theta )} =\P [M_\theta =x]\]
so \(X|_{\{\Theta =\theta \}}\eqd M_\theta \).
In the context of Definition 8.1.1 the ABC algorithm for generating samples from \(\Theta |_{\{X=x\}}\) is the following:
-
1. Sample \(\theta _0\) from the discrete distribution \(\Theta \).
-
2. Sample \(x_0\) from the discrete distribution \(M_\theta \eqd X|_{\{\Theta =\theta \}}\).
-
3. Then:
-
- if \(x\neq x_0\), go back to step one;
-
- if \(x=x_0\),accept \(\theta _0\) as a sample of \(\Theta |_{\{X=x\}}\).
-
This algorithm is precisely the strategy of our proof for Lemma 1.4.1, written as an algorithm and adapted to the special case of Definition 8.1.1. It generates a single sample of the distribution \(\Theta |_{\{X=x\}}\). We can run the algorithm again to obtain more samples.
The ABC algorithm outlined above requires only that we have the ability to take samples from discrete distributions with known probability mass functions. To handle cases with continuous priors and/or data, we also need to be able to sample from continuous distributions with known probability density functions. The modifications are as follows:
-
• If \(\Theta \) is continuous and \((M_\theta )\) is discrete then we can adopt the same algorithm, with the modification that in step 1 we must now sample from a continuous distribution rather than a discrete distribution.
-
• If \((M_\theta )\) is continuous then in step 3 we will have \(\P [x=x_0]=0\). In this case the simplest strategy is to fix some \(\eps >0\) and accept \(\theta _0\) as an approximate sample of \(\Theta |_{\{X=x\}}\) if \(|x-x_0|\leq \eps \).
This idea is based on (1.9), which stated that if \(\Theta |_{\{X=x\}}\) was to be defined then it should be defined to be the limit as \(\eps \to 0\) of \(\Theta |_{\{|X-x|\leq \eps \}}\). The terminology ‘Approximate’ Bayesian Computation comes from this step.
More complex strategies for comparing \(x\) and \(x_0\) can also be used, with the aim of focusing on the aspects of the data that are most important to us.
The ABC algorithm as described above has a serious drawback. In discrete cases the probability that \(x=x_0\) (in step 3) can be extremely small, meaning that we have to go around the loop \(1\to 2\to 3\to 1\) many times before we find a sample we can accept. In continuous cases, choosing \(\eps \) close to \(0\) obtains results in good approximation but introduces the same problem; that the probability of accepting the same \(x_0\) becomes small. To get handle this difficulty, various ways of sampling \(x_0\) (in step 2) have been developed to increase the acceptance probability, without changing the distribution of \(x_0\). One such method is ABC-MCMC, which uses ideas from Section 8.3 to sample \(x_0\). Another method is sequential ABC, where \(x_0\) is sampled as a perturbation in a carefully chosen direction from the (rejected) \(x_0\) in the previous iteration of the loop. We will not detail such methods here. They are very popular in some applications.