$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.