# MCEM algorithm in normal distribution

Consider $z_1,\ldots,z_n$ as a sample of observations of $Z$ and $y_1,\ldots,y_n$ the missing data,

where $Z\sim N(\mu,\sigma^2+1)$ and $Y\sim N(0,1)$.

i)Find the expression of $u_{p+1}$ and $\sigma^2_{p+1}$ given the current values $\mu_p$ and $\sigma^2_p$.

ii)Implement in R the EM algorithm for estimate $\mu$ and $\sigma^2$.Use as data a sample of $rnorm$

First I found that $Y\mid Z=z\sim N(\frac{z-\mu}{\sigma^2+1},\frac{\sigma^2}{\sigma^2+1})$. Then I found the likelihood

$$L\propto (\frac{\sigma^2+1}{\sigma^2})^{\frac{n}{2}}e^{-\frac{1}{2}\frac{\sigma^2+1}{\sigma^2}\sum(y-(\frac{z-\mu}{\sigma^2+1}))^2}$$ $$\log L=\frac{n}{2}\log(\sigma^2+1)-\frac{n}{2}\log(\sigma^2)-\frac{1}{2}\frac{\sigma^2+1}{\sigma^2}\sum (y-(\frac{z-\mu}{\sigma^2+1}))^2$$

So taking the derivative $$\frac{\partial}{\partial\mu}\log L=-\frac{\sigma^2+1}{\sigma^2}\frac{1}{\sigma^2+1}\sum (y-\frac{z-\mu}{\sigma^2+1})=-\frac{1}{\sigma^2}(\sum y-\frac{nz}{\sigma^2+1}+\frac{n\mu}{\sigma^2+1})$$ developing this I get $$\hat{\mu}=-(\sigma^2+1)\overline{Y}+z$$

This will be my $\mu_{p+1}$? I am studying through this book Casella.

THE EM Algorithm

1. Compute $Q(\theta\mid\hat{\theta}_{(m)},x)=E_{\hat{\theta}_{(m)}}[\log L^c(\theta\mid x,z)]$ where the expectation is with respect to $k(z\mid\hat{\theta}_m,x)=\frac{f(x,z\mid\theta)}{g(x\mid\theta)}$

2. Maximize $Q(\theta\mid\hat{\theta}_{(m)},x)$ in $\theta$ and take $\theta_{(m+1)}=\operatorname{argmax}\limits_\theta$$Q(\theta\mid \hat{\theta}_{(m)},x)$

EDIT: I think perhaps a mistake as I did, because if $\hat{\mu}$ is the $u_{p+1}$ it will not depending on the value of the previous iteration

EDIT2: I think there is no way to make the M step without using simulation, so confused with this topic, this would be a case of mixture of distributions, since Y and Z have different distributions?

• I did some caculation, I think $\mu_{p+1}=\mu_{p}=\mu$ for any iteration for this case. I don't know how to solve the variance. Hope somebody else can solve. Jun 19, 2015 at 0:50
• You shouldln't write $Z$~$N$; rather the proper usage is $Z\sim N$, coded as Z\sim N. That gets you proper spacing and no mismatch in fonts. Likewise \mid afford proper spacing in things like $\Pr(A\mid B)$ (as opposed to $\Pr(A|B)$, which is not correct usage). I also changed $log L$ to $\log L$, coded as \log L. That not only prevents italicization but also affords proper spacing in expressions like $a\log b$. I also cleaned up your "argmax". I await approval of the edits. ${}\qquad{}$ Jun 19, 2015 at 20:35
• (But my edits have so far not address sizes of parentheses, e.g. $\displaystyle\left(\sum_{i=1}^n\frac{A_i}{B_i}\right)$ as opposed to $\displaystyle(\sum_{i=1}^n\frac{A_i}{B_i})$. Just using \left( and \right) takes care of that. I don't use \left( and \right) in things like f(x) since there's no size issue then and it makes the code cluttered.) ${}\qquad{}$ Jun 19, 2015 at 20:37
• Also: Why write "Where $Z\sim N(\mu,\sigma^2+1)$ and $Y\sim N(0,1)$." with a capital initial "W" as if it were the beginning of a new sentence?? I've seen that only on the internet. I don't understand the appeal of such an absurdity. ${}\qquad{}$ Jun 19, 2015 at 20:40
• @askazy : I'm quite rusty on EM algorithms, but I'll take a look. Jun 19, 2015 at 20:58

Ok, let me show how I did it for $\mu$, but if i am wrong, i hope someone can point out.

First let us compute the complete likelyhood (i.e the likelyhood for observed and missing together, or some people call it as augmented likelyhood).This is not difficult.

$L^c(\mu|z,y)\propto exp\left \{ -\frac{1}{2}[\sum_{i=1}^{n}(z_i-\mu)^2 +\sum_{i=1}^ny_{i}^2]\right \}$

$z_i$ is observed value with $N(\mu,\sigma^2+1)$ distribution, here we suppose $\sigma$ is known. $y_i$ are missing values. We omitted all constants for the likelihood function.

Next we calculate the likelihood for the observed values.

$L(\mu|z)\propto exp\left \{ -\frac{1}{2}\sum_{i=1}^{n}(z_i-\mu)^2 \right \}$

The condition pdf of missing $y_i$ conditional on observed $z_i$ and $\mu$ is calculated by $L^c(\mu|z,Y)/L(\mu|z)$ which is: $exp(-\frac{1}{2}\sum_{i}^ny_i^2$). This is the product of i.i.d distribution of $N(0,1)$

It seems like we even don't need this step to show that missing values have a $N(0,1)$ distribution,

Anyway, now we know that conditional on observed $z_i$ and $\mu$, missing value $y_i$ has a $N(0,1)$ distribution.

Next step, we need to calculate the expected value of log likelihood of complete (or augmented) likelihood function conditional on an initially guessed $\mu_0$ and observed values $z_i$. i.e $E[logL^c(\mu|z,Y)|\mu_0,z]$

Further we can write this expectation as

$E\left \{-\frac{1}{2}[\sum_{i=1}^{n}(z_i-\mu)^2 +\sum_{i=1}^ny_{i}^2]|\mu_0,z\right \}$

Since $Z_i$ are all observed value. The expected value for $-\frac{1}{2}\sum_{i=1}^{n}(z_i-\mu)^2$ is just itself. i.e we treat the first part as constant.

So, for the second part we plug in an initial guess $\mu_0$, then we can use a condition we get before, i.e conditional on $\mu$ (here is $\mu_0$) and observed $z_i$, $Y$ has a $N(0,1)$ distribution.

So, $E \left \{ -\frac{1}{2}[\sum_{i=1}^ny_i^2]|\mu_0,z\right \}=E \left \{ -\frac{1}{2}[\sum_{i=1}^n(y_i-\mu_0+\mu_0)^2]|\mu_0,z\right \}=E \left \{ -\frac{1}{2}[\sum_{i=1}^n(y_i-\mu_0)^2+2\sum_{i=1}^n(y_i-\mu_0)\mu_0+\mu_0^2]|\mu_0,z\right \}=-\frac{1}{2}(n+\mu_0^2)$

Note $E(y_i-\mu_0)^2$ is the variance of $y_i$ and $\sum_{i=1}^n(y_i-\mu_0)=0$

So the Expectation for the complete log likelihood

$E\left \{-\frac{1}{2}[\sum_{i=1}^{n}(z_i-\mu)^2 +\sum_{i=1}^ny_{i}^2]|\mu_0,z\right \}=-\frac{1}{2}[\sum_{i=1}^{n}(z_i-\mu)^2+n+\mu_0^2]$

Next take partial derivative of the expectation in term of $\mu$ and set it to zero

we get $\mu=\bar{z}$, we even don't need a initial $\mu_0$ for this case.

Anyway, I think I have not totally understand the EM algorithm yet,the solution might be totally wrong.