MCEM algorithm in normal distribution

Question

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?

Answer 1

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.