Conditional distribution for Gibbs sampling for Gaussian mixture

Question

If we draw $n$ i.i.d. points $x_1,x_2,\dots,x_n$ from the following Gaussian mixture: $$ \frac 12 \mathcal N(x \mid \mu_1,1) + \frac 12 \mathcal N(x\mid \mu_2,1) $$ and the prior $p(\mu_1 , \mu_2 )$ is $\mathcal N (\mu_1 \mid 0, \tau )\mathcal N (\mu_2 |0, \tau )$, where $\tau$ is a fixed constant.

What would be the conditional distributions from which we can sample in the Gibbs Sampling Algorithm? How do we sample from this kind of distribution?

Answer 1

The Gibbs steps for a mixture model are to be found in all papers and books addressing Bayesian inference on mixtures, from our early paper with Diebolt (1990) to the reference book of Sylvia Frühwirth-Schnatter. See for instance our own review with Jean-Michel Marin and Kerrie Mengersen. The following is taken verbatim from our book Bayesian Essentials with R and works with the marginal posterior on the component allocation auxiliary vector ${\mathbf Z}$:

With each $x_i$ is associated a missing variable $z_i$ that indicates "its" component, i.e. the index $z_i$ of the distribution from which it was generated. Formally, this means that we have a hierarchical structure associated with the model: $$ z_i|{\mathbf p}\sim\mathscr{M}_k(p_1,\ldots, p_k) $$ and $$ x_i|z_i,{\mathbf \theta} \sim f(\cdot|\theta_{z_i})\,. $$ In the case of the unknown mean Normal mixture, if we take two different independent Normal priors on both means, $$ \mu_1 \sim \mathscr{N}(0,4)\,,\quad\mu_2\sim \mathscr{N}(2,4)\,, $$ the posterior weight of a given allocation vector ${\mathbf z}$ is \begin{align*} \omega\left({\mathbf z}\right) \propto &\sqrt{(n_1+1/4)(n-n_1+1/4)}\,p^{n_1} (n_1-p)^{n-l} \,\\ &\times\exp\left\{-[(n_1+1/4)\hat s_1\left({\mathbf z}\right) + n_1\{\bar x_1\left({\mathbf z}\right)\}^2/4]/2 \right\} \\ &\times\exp\left\{-[(n-n_1+1/4)\hat s_2\left({\mathbf z}\right) + (n-n_1)\{\bar x_2\left({\mathbf z}\right)-2\}^2/4]/2 \right\}, \end{align*} $$\displaystyle \bar x_1\left({\mathbf z}\right)=\frac{1}{n_1}\sum_{i=1}^n\mathbb{I}_{z_i=1}x_i,\quad \bar x_2\left({\mathbf z}\right)=\frac{1}{n-n_1}\sum_{i=1}^n\mathbb{I}_{z_i=2}x_i\,,$$ $$\displaystyle \hat s_1\left({\mathbf z}\right)=\sum_{i=1}^n\mathbb{I}_{z_i=1}\left(x_i-\bar x_1\left({\mathbf z}\right)\right)^2,\quad \hat s_2\left({\mathbf z}\right)=\sum_{i=1}^n\mathbb{I}_{z_i=2}\left(x_i-\bar x_2\left({\mathbf z}\right)\right)^2$$ (if we set $\bar x_1\left({\mathbf z}\right)=0$ when $n_1=0$ and $\bar x_2\left({\mathbf z}\right)=0$ when $n-n_1=0$). Implementing this derivation in R is quite straightforward:

  omega=function(z,x,p){
  n=length(x)
  n1=sum(z==1);n2=n-n1
  if (n1==0) xbar1=0 else xbar1=sum((z==1)*x)/n1
  if (n2==0) xbar2=0 else xbar2=sum((z==2)*x)/n2
  ss1=sum((z==1)*(x-xbar1)^2)
  ss2=sum((z==2)*(x-xbar2)^2)
  return(sqrt((n1+.25)*(n2+.25))*p^n1*(1-p)^n2*
     exp(-((n1+.25)*ss1+(n2+.25)*ss2)/2)*
     exp(-(n1*xbar1^2+n2*xbar2)/8))
  }

leading for instance to

> omega(z=sample(1:2,4,rep=TRUE),x=plotmix(n=4,plot=FALSE)$samp,p=.8)
[1] 0.0001781843
> omega(z=sample(1:2,4,rep=TRUE),x=plotmix(n=4,plot=FALSE)$sample,p=.8)
[1] 5.152284e-09

Note that the omega function is not and cannot be normalized, so the values must be interpreted on a relative scale.

Here is a further excerpt from Monte Carlo Statistical Methods (2004, p.342) that spells out the Gibbs sampler in full details:

Consider a normal mixture with two components with equal known variance and fixed weights, $$ p\,\mathcal{N}(\mu_1,\sigma^2) + (1-p)\,\mathcal{N}(\mu_2,\sigma^2) \,. $$ We assume in addition a normal $\mathcal{N}(0,10\sigma^2)$ prior distribution on both means $\mu_1$ and $\mu_2$. Generating directly from the posterior associated with a sample $\mathbf{x} = (x_1,\ldots,x_n)$ from this normal mixture quickly turns impossible, as discussed for instance in Diebolt and Robert (1994) and Celeux, Hurn and Robert (2000), because of a combinatoric explosion in the number of calculations, which grow as $\mathcal{O}(2^n)$.

As for the EM algorithm, a natural completion of $(\mu_1,\mu_2)$ is to introduce the (unobserved) component indicators $z_i$ of the observations $x_i$, namely, $$ \mathbb{P}(Z_i=1) = 1-\mathbb{P}(Z_i=2) = p \qquad\mbox{and}\qquad X_i|Z_i=k\sim\mathcal{N}(\mu_k,\sigma^2) \,. $$ The completed distribution is thus \begin{align*} \pi(\mu_1,\mu_2,\mathbf{z}|\mathbf{x}) &\propto \exp\{-(\mu_1^2+\mu_2^2)/20 \sigma^2 \}\, \prod_{z_i=1} p \exp\{-(x_i-\mu_1)^2/2\sigma^2 \}\,\times\\ &\prod_{z_i=2} (1-p) \exp\{-(x_i-\mu_2)^2/2\sigma^2 \}\,. \end{align*} Since $\mu_1$ and $\mu_2$ are independent, given $(\mathbf{z},\mathbf{x})$, with distributions $(j=1,2)$, the conditional distributions are $$\mathcal{N} \left( \sum_{z_i=j} x_i \big/ \left( .1 + n_j\right), \sigma^2 \big/ \left( .1 + n_j\right) \right)\,,$$ where $n_j$ denotes the number of $z_i$'s equal to $j$. Similarly, the conditional distribution of $\mathbf{z}$ given $(\mu_1,\mu_2)$ is a product of binomials, with $$ {\mathbb P}(Z_i=1|x_i,\mu_1,\mu_2) = \frac{p \exp\{-(x_i-\mu_1)^2/2\sigma^2 \} }{ p \exp\{-(x_i-\mu_1)^2/2\sigma^2 \} + (1-p) \exp\{-(x_i-\mu_2)^2/2\sigma^2 \}} \,. $$ Figure above illustrates the behavior of the Gibbs sampler in that setting, with a simulated dataset of $500$ points from the $.7 \mathcal{N}(0,1) +.3 \mathcal{N}(2.7,1)$ distribution. The representation of the MCMC sample after $15,000$ iterations is quite in agreement with the posterior surface, represented via a grid on the $(\mu_1,\mu_2)$ space and some contours; while it may appear to be too concentrated around one mode, the second mode represented on this graph is much lower since there is a difference of at least $50$ in log-posterior values. However, the Gibbs sampler may also fail to converge, as described in Diebolt and Robert (1994) and illustrated in Figure above. When initialized at a local mode of the likelihood, the magnitude of the moves around this mode may be too limited to allow for exploration of further modes (in a reasonable number of iterations).