Solution (source code)

= Solution

Introduce a categorical <latent variable> $Z_i$ for each observation, with $P(Z_i=j)=\alpha_j$, and let $z_{ij}=\mathbf1_{\{Z_i=j\}}$. Conditional on $Z_i=j$, $X_i$ has a <normal distribution> of mean $\mu_j$ and the shared variance $\sigma^2$. The complete-data <log-likelihood function> is
$$
\ell_c(\theta;x,z)=\sum_{i=1}^N\sum_{j=1}^kz_{ij}\log\alpha_j-\frac N2\log(2\pi\sigma^2)-\frac1{2\sigma^2}\sum_{i=1}^N\sum_{j=1}^kz_{ij}(x_i-\mu_j)^2.
$$
Initialize finite means, positive weights summing to one, and a positive common variance. At iteration $r$, the <Bayes' theorem> gives the E-step <mixture responsibilities>
$$
\boxed{r_{ij}=P(Z_i=j\mid x_i,\theta^{(r)})=\frac{\alpha_j^{(r)}\exp[-(x_i-\mu_j^{(r)})^2/(2\sigma^{2(r)})]}{\sum_{\ell=1}^k\alpha_\ell^{(r)}\exp[-(x_i-\mu_\ell^{(r)})^2/(2\sigma^{2(r)})]}.}
$$
The common normal-density factor cancels. The observations' conditional labels are independent, and replacing $z_{ij}$ by their conditional expectations $r_{ij}$ gives
$$
Q(\theta\mid\theta^{(r)})=\sum_{i,j}r_{ij}\log\alpha_j-\frac N2\log\sigma^2-\frac1{2\sigma^2}\sum_{i,j}r_{ij}(x_i-\mu_j)^2+\text{constant}.
$$
Set $n_j=\sum_i r_{ij}$; because $\sum_jr_{ij}=1$, $\sum_jn_j=N$. The constrained weight maximization uses a <Lagrange multiplier> $\eta$: differentiating $\sum_jn_j\log\alpha_j-\eta(\sum_j\alpha_j-1)$ gives $\alpha_j=n_j/\eta$, and summing fixes $\eta=N$. Differentiating the residual quadratic in each mean gives $\sum_i r_{ij}(x_i-\mu_j)=0$. Hence the first two M-step updates are
$$
\boxed{\alpha_j^{(r+1)}=\frac{n_j}{N},\qquad\mu_j^{(r+1)}=\frac{\sum_i r_{ij}x_i}{n_j}.}
$$
With these new means define $R=\sum_{i,j}r_{ij}(x_i-\mu_j^{(r+1)})^2$. The remaining objective in $u=\sigma^2$ is $-N\log u/2-R/(2u)$, whose derivative is $-N/(2u)+R/(2u^2)$. For $R>0$ it increases up to $u=R/N$ and decreases afterwards. Thus the shared-variance update is
$$
\boxed{\sigma^{2(r+1)}=\frac1N\sum_{i=1}^N\sum_{j=1}^kr_{ij}(x_i-\mu_j^{(r+1)})^2.}
$$
The denominator is $N$, not $N-k$ and not a separate count for each component: this is <maximum likelihood estimation> with one common variance. These are the <EM for Gaussian mixtures with a common variance> updates. Recompute responsibilities with the new parameters and repeat until the observed-data likelihood changes sufficiently little.

The monotonicity claim also follows directly. For positive old weights, the likelihood ratio for observation $i$ is
$$
\frac{\sum_j\alpha_j'f(x_i;\mu_j',\sigma'^2)}{\sum_j\alpha_j^{(r)}f(x_i;\mu_j^{(r)},\sigma^{2(r)})}=\sum_jr_{ij}\frac{\alpha_j'f(x_i;\mu_j',\sigma'^2)}{\alpha_j^{(r)}f(x_i;\mu_j^{(r)},\sigma^{2(r)})}.
$$
Apply the <Jensen inequality> to the logarithm and sum over observations. This yields
$$
\ell(\theta')-\ell(\theta^{(r)})\geq Q(\theta'\mid\theta^{(r)})-Q(\theta^{(r)}\mid\theta^{(r)})\geq0
$$
for the maximizing M step. This proves likelihood increase, not convergence to the global maximum. Use several initializations to reduce sensitivity to local maxima. Compute responsibilities with a <log-sum-exp function> normalization to avoid numerical underflow.

If $n_j=0$, that component's updated weight is zero and its mean is unidentified; retain its old mean or remove the empty component. Strictly positive interior initialization avoids this case in exact arithmetic while the variance stays positive. If $R=0$, there is no positive finite variance maximizing this M step: the objective tends to infinity as $\sigma^2\downarrow0$. Such boundary degeneracy needs a variance constraint or regularization, as discussed in the next solution. For $k=1$ the formulas reduce to the ordinary sample mean and the maximum-likelihood variance $N^{-1}\sum_i(x_i-\overline x)^2$.