# Conditional distribution for Gibbs sampling for Gaussian mixture

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?

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))
}


> 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).

• where and how have you used the mean values (0 and 2) from the prior of means? how do you define this conditional distribution? – NerdHardAtWork Nov 9 '18 at 1:06
• if we were doing the same thing for Metropolis-Hastings, how would we choose the proposal distribution? – NerdHardAtWork Nov 9 '18 at 1:07
• (1) In this resolution, both means are integrated out, this is the marginal distribution of the vector $z$; (2) Metropolis-Hastings starts from $\omega({\mathbf z})$ as a target and makes a proposal for ${\mathbf z}$. – Xi'an Nov 9 '18 at 4:36
• Thanks. Is there any R code for the Gibbs sampler in the second part? – user261225 Nov 27 '20 at 16:08