EM with Binomial $X_1, X_2 \dots X_n$ represent independently observed Bernoulli random variables.
$Z_1, Z_2, \dots Z_n$ are unobserved.
$Z_i | \theta_i \sim N(\theta_i,1)$
$\theta_i \sim N(\epsilon, \sigma^2)$
$X_i = \begin{cases} 0 &  Z_i \leq u \\ 1 & Z_i > u \end{cases}$ where $u$ is known.
Show that the EM sequence for estimation of $\epsilon$ is given by:
$\epsilon^{t+1} = \frac{1}{n} \sum_{i=1}^n E[Z_i|X_i, \epsilon^t, \sigma^2]$
$\Phi$ represents CDF of standard normal.
Attempt
$Z_i \sim N(\epsilon, \sigma^2+1)$ 
\begin{align*}
f(X,Z| \epsilon, \sigma^2) &= \prod_{i=1}^n p_i^{x_i} (1-p_i)^{1-x_i}\\
p_i &= P(Z_i \leq u)\\
&=\Phi(\frac{u-\epsilon}{\sqrt{\sigma^2+1}})
\end{align*}
Complete log likeihood = $\sum_{i=1}^nx_i \log p_i + (n-\sum_i x_i) \log(1-p_i)$
I am not able to proceed further.
 A: Define $p = P(Z_i > u) = P( \frac{Z-\mathcal{\epsilon}}{\sqrt{\sigma^2+1}} > \frac{u-\mathcal{\epsilon_i}}{\sqrt{\sigma^2+1}}) = \phi(\frac{\mathcal{\epsilon_i}-u}{\sigma{\sigma^2+1}})$
The incomplete data log likelihood:
$$
\begin{align*}
L(\theta | X) &= \prod_{i=1}^n p^{x_i}(1-p)^{1-x_i}\\
&= p^{\sum_i x_i}(1-p)^{n-\sum_i x_i}
\end{align*}
$$
Complete data log likehood:
$$
\begin{align*}
f(Z|\mathcal{\epsilon}, \sigma^2) &= (2\pi(\sigma^2+1))^{-\frac{n}{2}} \exp(-\sum_i\frac{(z_i-\mathcal{\epsilon})^2}{2(\sigma^2+1)})\\
\log L(\mathcal{\epsilon} | Z,\sigma^2) & = -\frac{n}{2} \log(2\pi(\sigma^2+1)) -\sum_i\frac{(z_i-\mathcal{\epsilon})^2}{2(\sigma^2+1)}\\
E[\log L(\mathcal{\epsilon} | Z,\sigma^2) ] &= -\frac{n}{2} \log(2\pi(\sigma^2+1)) -\frac{1}{2(\sigma^2+1)} E[\sum_i z_i^2 -2\mathcal{\epsilon}\sum_i z_i + n\mathcal{\epsilon^2}]
\end{align*}
$$
where the expectation is with repsect to the conditonal distribution of $Z$ having observed $X$, $k(z|x,\theta) = h(x,z|\theta)/g(x | \theta)$
$$
\begin{align*}
Q(\mathcal{\epsilon}| \mathcal{\epsilon_0}, \sigma^2, X) &= -\frac{n}{2} \log(2\pi(\sigma^2+1)) -\frac{1}{2(\sigma^2+1)} E[\sum_i z_i^2 -2\mathcal{\epsilon}\sum_i z_i + n\mathcal{\epsilon^2}]\\
\frac{\partial Q}{\partial \mathcal{\epsilon}} &= -\frac{1}{2(\sigma^2+1)} \big( -2 \sum_iE[ z_i | x_i] + 2n\mathcal{\epsilon} \big) = 0\\
\implies \mathcal{\epsilon^{(i)}} = \frac{1}{n} E[z_i|x_i,\mathcal{\epsilon^{i-1},\sigma^2}]
\end{align*}
$$
