# Estimation of parameter when having two sources of errors

First of all sorry for the bad title. I picture this is a standard problem, i just don't know how to deal with it

I have the following model

\begin{equation} Y_{ij} = \mu + E_i + \varepsilon_{ij} \end{equation} Where $1 \leq i \leq N$ and $1 \leq j \leq M$. And $E_i \sim \mathcal{N}(0,\sigma^2_E$) and $\varepsilon_{ij}\sim \mathcal{N}(0,\sigma^2_\varepsilon)$ are error terms all independent among each other. Note that $E_i$ appear in many equations and $\varepsilon_{ij}$ in only one.

I would like to know how is the standard way to estimate $\mu$ and the variance of the estimation, after observing all $Y_{ij}$.

My first approach is to estimate by the mean \begin{equation} \hat{\mu} =\frac{1}{MN} \sum_{i=1}^{N}\sum_{j=1}^M Y_{ij} \end{equation} I am able to calculate the variance of $\hat{\mu}$, that yields \begin{equation} \text{Var}(\hat{\mu}) = \frac{\sigma^2_{\varepsilon}}{MN} + \frac{\sigma^2_E}{N} \end{equation} I can manage to estimate $\sigma^2_\varepsilon$ with \begin{equation} \hat{\sigma^2_\varepsilon} = \frac{1}{N(M-1)}\sum_{i=1}^M\sum_{j=1}^N(Y_{ij} - \bar{Y_i})^2 \end{equation} But I am not sure how to estimate $\sigma^2_E$.