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?


1 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:

  if (n1==0) xbar1=0 else xbar1=sum((z==1)*x)/n1
  if (n2==0) xbar2=0 else xbar2=sum((z==2)*x)/n2

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 \}} \,. $$ enter image description here 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. enter image description here 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).

  • $\begingroup$ where and how have you used the mean values (0 and 2) from the prior of means? how do you define this conditional distribution? $\endgroup$ Nov 9, 2018 at 1:06
  • $\begingroup$ if we were doing the same thing for Metropolis-Hastings, how would we choose the proposal distribution? $\endgroup$ Nov 9, 2018 at 1:07
  • $\begingroup$ (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}$. $\endgroup$
    – Xi'an
    Nov 9, 2018 at 4:36
  • $\begingroup$ Thanks. Is there any R code for the Gibbs sampler in the second part? $\endgroup$
    – user261225
    Nov 27, 2020 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.