How to derive Gibbs sampling?

Question

I'm actually hesitating to ask this, because I'm afraid I will be referred to other questions or Wikipedia on Gibbs sampling, but I don't have the feeling that they describe what's at hand.

Given a conditional probability $p(x|y)$: $$ \begin{array}{c|c|c} p(x|y) & y = y_0 & y = y_1 \\ \hline x = x_0 & \tfrac{1}{4} & \tfrac{2}{6} \\ \hline x = x_1 & \tfrac{3}{4} & \tfrac{4}{6} \\ \end{array} $$

And a conditional probability $p(y|x)$: $$ \begin{array}{c|c|c} p(y|x) & y = y_0 & y = y_1 \\ \hline x = x_0 & \tfrac{1}{3} & \tfrac{2}{3} \\ \hline x = x_1 & \tfrac{3}{7} & \tfrac{4}{7} \\ \end{array} $$

We can uniquely come up with the joint probability $f_{unique}=p(x,y)$:

$$ \begin{array}{c|c|c|c} p(x,y) & y = y_0 & y = y_1 & p(x) \\ \hline x = x_0 & a_0 & a_1 & c_0 \\ \hline x = x_1 & a_2 & a_3 & c_1 \\ \hline p(y) & b_0 & b_1 & \\ \end{array} $$

Because, although we have $8$ unknowns, we have more ($4*2+3$) linear equations:

$ a_0+a_1+a_2+a_3=1 \\ b_0+b_1 = 1 \\ c_0+c_1 = 1 $

As well as:

$ \tfrac{1}{4} b_0 = a_0 \\ \tfrac{3}{4} b_0 = a_2 \\ \tfrac{2}{6} (1-b_0) = a_1 \\ \tfrac{4}{6} (1-b_0) = a_3 \\ \tfrac{1}{3} c_0 = a_0 \\ \tfrac{2}{3} c_0 = a_1 \\ \tfrac{3}{7} (1-c_0) = a_2 \\ \tfrac{4}{7} (1-c_0) = a_3 $

It's quickly solved by $c_0=\tfrac{3}{4}b_0$, $\tfrac{2}{3}c_0=a_1$. Namely by equating $\tfrac{2}{4}b_0=a_1$ with $\tfrac{2}{6}(1-b_0)=a_1$. This gives $b_0=\tfrac{2}{5}$ and the rest follows.

$$ \begin{array}{c|c|c|c} p(x,y) & y = y_0 & y = y_1 & p(x) \\ \hline x = x_0 & \tfrac{1}{10} & \tfrac{2}{10} & \tfrac{3}{10} \\ \hline x = x_1 & \tfrac{3}{10} & \tfrac{4}{10} & \tfrac{7}{10} \\ \hline p(y) & \tfrac{4}{10} & \tfrac{6}{10} & \\ \end{array} $$

So, now we go to the continuous case. It is imaginable to go to intervals and keep the above structure in-tact (with more equations than unknowns). However, what happens when we go to (point) instances of random variables? How does sampling

$$ x_a \sim p(x|y=y_b) \\ y_b \sim p(y|x=x_a) $$

iteratively, lead to $p(x,y)$? Equivalent to the constraint $a_0 + a_1 + a_2 + a_3=1$, how does it ensure $\int_X \int_Y p(x,y) dy dx = 1$ for example? Likewise with $\int_Y p(y|x)dy=1$. Can we write down the constraints and derive Gibbs sampling from first principles?

So, I'm not interested in how to perform Gibbs sampling, which is simple, but I'm interested in how to derive it, and preferably how to prove that it works (probably under certain conditions).

Answer 1

Computing a joint distribution from conditional distributions in general is very difficult. If the conditional distributions are chosen arbitrarily, a common joint distribution might not even exist. In this case, even showing that the conditional distributions are consistent is generally difficult. One result that might be used for deriving a joint distribution is Brook's lemma, $$ \frac{p(\mathbf{x})}{p(\mathbf{x}')} = \prod_i \frac{p(x_i \mid \mathbf{x}_{<i}, \mathbf{x}'_{>i})}{p(x_i' \mid \mathbf{x}_{<i}, \mathbf{x}'_{>i})},$$ by choosing a fixed state $\mathbf{x}'$, although I have never successfully used it myself for that purpose. For more on that topic, I would look at Julian Besag's work.

To prove that Gibbs sampling works, however, it's better to take a different route. If a Markov chain implemented by a sampling algorithm has distribution $p$ as invariant distribution, and is irreducible and aperiodic, then the Markov chain will converge to that distribution (Tierney, 1994).

Gibbs sampling will always leave the joint distribution invariant from which the conditional distributions were derived: Roughly, if $(x_0, y_0) \sim p(x_0, y_0)$ and we sample $x_1 \sim p(x_1 \mid y_0)$, then

$$(x_1, y_0) \sim \int p(x_0, y_0) p(x_1 \mid y_0) \, dx_0 = p(x_1 \mid y_0) p(y_0) = p(x_1, y_0).$$

That is, updating $x$ by conditionally sampling does not change the distribution of the sample.

However, Gibbs sampling is not always irreducible. While we can always apply it without breaking things (in the sense that if we already have a sample from the desired distribution it will not change the distribution), it depends on the joint distribution whether Gibbs sampling will actually converge to it (a simple sufficient condition for irreducibility is that the density is positive everywhere, $p(\mathbf{x}) > 0$).