1
$\begingroup$

I'm solving an exercise from Casella's book, Introduction to Monte Carlo with R and he askes us to build a Gibbs sampler for the following mixture of Normals

\begin{align} p\mathcal{N}(\mu_{1},\sigma^{2}) + (1-p)\mathcal{N}(\mu_{2},\sigma^{2}) \end{align} such that the prior for both means is $\mu_{j}\sim\mathcal{N}(0,\sigma^{2}v^{2})$. We proceed by doing data augmentation and defining the latent variable $Z_{i}\sim Bernoulli(p)$ so we can properly divide the posterior distribution and finally get to this simpler form:

\begin{align} f(\mu_{1},\mu_{2},\textbf{z}|\textbf{x}) &\propto \exp \left\{-\frac{(\mu_{1}^{2}+\mu_{2}^{2})}{2v^{2}\sigma^2}\right\} \prod_{i:z_{i}=1}p\exp\left\{-\frac{(x_{i}-\mu_{1})^2}{2\sigma^2}\right\} \prod_{i:z_{i}=0}(1-p)\exp\left\{-\frac{(x_{i}-\mu_{2})^2}{2\sigma^2}\right\} \end{align}

Casella askes us to verify that the full conditionals of each $\mu_{j}$ are $\mathcal{N}\left(\frac{v^2\sum_{i:z_{i}=j} x_i}{n_{j}v^2+1},\frac{v^2\sigma^2}{n_{j}v^2+1}\right)$ such that $n_{j}$ is simply the number of observations that are in that "class". I rewrote the distribution as

\begin{align} f(\mu_{1},\mu_{2},\textbf{z}|\textbf{x}) & \propto p^{n_{1}}\exp\left\{-\frac{\mu_{1}^2+v^2\sum_{i:z_{i}=0}(x_{i}-\mu_{1})^2}{2v^2\sigma^2}\right\} (1-p)^{n_{2}}\exp\left\{-\frac{\mu_{2}^{2}+v^2\sum_{i:z_{i}=1}(x_{i}-\mu_{2})^2}{2v^2\sigma^2}\right\} \end{align}

Looking at the exponential part especifically \begin{align} \exp\left\{-\frac{\mu_{1}^2+v^2\sum_{i:z_{i}=0}(x_{i}-\mu_{1})^2}{2\sigma^2}\right\} &=\exp\left\{-\frac{\mu_{1}^2+v^2\sum_{i:z_{i}=0}x_{i}^2-2v^2\mu_{1}\sum_{i:z_{i}=0}x_{i} +v^2n_{1}\mu_{1}^2}{2\sigma^2}\right\} \\ &=\exp\left\{-\frac{(1+v^2n_{1})\mu_{1}^2-2v^2\mu_{1}\sum_{i:z_{i}=0}x_{i}+v^2\sum_{i:z_{i}=0}x_{i}^2}{2v^2\sigma^2}\right\} \\ &=\exp\left\{-\frac{1+v^2n_{1}}{2v^2\sigma^2}\left(\mu_{1}^2-\frac{2v^2\mu_{1}\sum_{i:z_{i}=0}x_{i}}{1+v^2n_{1}}+\frac{v^2\sum_{i:z_{i}=0}x_{i}^2}{1+v^2n_{1}}\right)\right\} \end{align} Completing the square we get \begin{align} \exp\left\{-\frac{1+v^2n_{1}}{2v^2\sigma^2}\left[\left(\mu_{1}-\frac{v^2\sum_{i:z_{i}=j} x_i}{n_{j}v^2+1}\right)^2 +\frac{v^2\sum_{i:z_{i}=0}x_{i}^2}{1+v^2n_{1}}-\left(\frac{v^2\sum_{i:z_{i}=j} x_i}{n_{j}v^2+1}\right)^2\right] \right\} \end{align}

This would give me the mean and variance I want if it wasn't for the final terms. As for conditioning on $\textbf{Z}$, it would simply cancel the $p$ powers. What am i doing wrong? I would appreciate an insight.

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.