Gibbs sampling

In mathematics and physics, Gibbs sampling is an algorithm to generate a sequence of samples from the joint probability distribution of two or more random variables. The purpose of such a sequence is to approximate the joint distribution (i.e. to generate a histogram), or to compute an integral (such as an expected value). Gibbs sampling is a special case of the Metropolis-Hastings algorithm, and thus an example of a Markov chain Monte Carlo algorithm. The algorithm is named after the physicist J. W. Gibbs, in reference to an analogy between the sampling algorithm and statistical physics. The algorithm was devised by Geman and Geman (citation below), some eight decades after the passing of Gibbs, and is also called the Gibbs sampler.

Gibbs sampling is applicable when the joint distribution is not known explicitly, but the conditional distribution of each variable is known. The Gibbs sampling algorithm is to generate an instance from the distribution of each variable in turn, conditional on the current values of the other variables. It can be shown (see, for example, Gelman et al.) that the sequence of samples comprises a Markov chain, and the stationary distribution of that Markov chain is just the sought-after joint distribution.

Gibbs sampling is particularly well-adapted to sampling the posterior distribution of a Bayesian network, since Bayesian networks are typically specified as a collection of conditional distributions. BUGS (link below) is a program for carrying out Gibbs sampling on Bayesian networks.

Background
Gibbs sampling is similar to Metropolis-Hastings sampling but the value is always accepted ($$\left.\alpha = 1\right.$$). The point of Gibbs sampling is that given a bivariate distribution it is simpler to sample from a conditional distribution rather than integrating over a joint distribution. Suppose we want to sample $$\left.k\right.$$ values of $$\left.x\right.$$ from a joint distribution $$\left.p(x, y)\right.$$. We begin with a value of $$\left.y_0\right.$$ and sample $$\left.x\right.$$ by $$x_i \sim p\left(x | y = y_{i-1}\right)$$. Once that value of $$\left.x\right.$$ is calculated, repeat by sampling for the next $$\left.y\right.$$: $$y_i \sim p\left(y | x = x_i\right)$$.

Implementation
Suppose that a sample $$\left.X\right.$$ is taken from a distribution depending on a parameter vector $$\theta \in \Theta \,\!$$ of length $$\left.d\right.$$, with prior distribution $$g(\theta_1, \ldots, \theta_d)$$. It may be that $$\left.d\right.$$ is very large and that numerical integration to find the marginal densities of the $$\left.\theta_i\right.$$ would be computationally expensive. Then an alternative method of calculating the marginal densities is to create a Markov chain on the space $$\left.\Theta\right.$$ by repeating these two steps:


 * 1) Pick a random index $$1 \leq j \leq d$$
 * 2) Pick a new value for $$\left.\theta_j\right.$$ according to $$g(\theta_1, \ldots, \theta_{j-1} , \, \cdot \, , \theta_{j+1} , \ldots , \theta_d )$$

These steps define a reversible markov chain with the desired invariant distribution $$\left.g\right.$$. This can be proven as follows. Define $$x \sim_j y$$ if $$\left.x_i = y_i\right.$$ for all $$j \neq i$$ and let $$\left.p_{xy}\right.$$ denote the probability of a jump from $$x \in \Theta$$ to $$y \in \Theta$$. Then, for $$x \sim_j y$$ the transition probabilities are


 * $$p_{xy} = \frac{1}{d}\frac{g(y)}{\sum_{z \in \Theta: z \sim_j y} g(z) } $$

and $$\left.p_{xy} = 0\right.$$ otherwise. So



g(x) p_{xy} = \frac{1}{d}\frac{ g(x) g(y)}{\sum_{z \in \Theta: z \sim_j y} g(z) } = \frac{1}{d}\frac{ g(y) g(x)}{\sum_{z \in \Theta: z \sim_j x} g(z) } = g(y) p_{yx} $$

so the detailed balance equations are satisfied showing that the chain is reversible, and that it has invariant distribution $$\left.g\right.$$.

In practice, the suffix $$\left.j\right.$$ is not chosen at random, and the chain cycles through the suffixes in order. In general this gives a non-reversible chain, but it will still have the desired invariant distribution (as long as the chain can access all states under the fixed ordering).