Why did log-likelihood decrease in EM in this step-by-step example?

Question

I'm using Expectation Maximization algorithm to determine the parameters of Gaussian distributions in a mixture. To get a better understanding of the algorithm, I executed it manually step by step on a small example, following this procedure:

The problem is that in the second iteration log-likelihood decreased.

These are my data: $$X = \begin{bmatrix}1 \\ 2\\ 3\end{bmatrix}$$ The initial means, variances and mixing coefficients are set to:

$$\mu = \begin{bmatrix}1 & 3\end{bmatrix}\quad \pi = \begin{bmatrix}0.5 & 0.5\end{bmatrix}\quad \Sigma=\begin{bmatrix}1 & 1\end{bmatrix}$$

So, I first calculated the initial log-likelihood:

$$L_0 = -4.3893$$

Then I did the E step and calculated $\gamma(z_{nk})$, which I'll shorten as $\gamma_{n,k}$, like this $\gamma_{n,k}=\pi_kN(x_n|\mu_k,\Sigma_k)$, and got:

$$\Gamma = \begin{bmatrix}0.881 & 0.119 \\ 0.5 & 0.5 \\ 0.119 & 0.881\end{bmatrix}$$

Then, I peformed the M step: $$N_1 = 1.5\quad N_2= 1.5$$ $$\pi = \begin{bmatrix}0.5 & 0.5\end{bmatrix}\quad \mu^{new}=\begin{bmatrix}1.492 & 2.508\end{bmatrix}\quad \Sigma=\begin{bmatrix}0.409 & 0.409\end{bmatrix}$$ The new log-likelihood was $L_1=-3.6766$, which was greater than $L_0$.

Then, I did another iteration. In the E step I got: $$\Gamma=\begin{bmatrix}0.9979 & 0.0021 \\ 0.5 & 0.5 \\ 0.0021 & 0.9979 \end{bmatrix}$$ In the M step I got: $$N_1=1.5\quad N_2=1.5$$ $$\pi = \begin{bmatrix}0.5 & 0.5\end{bmatrix}\quad \mu^{new}=\begin{bmatrix}1.3361 & 2.6639\end{bmatrix}\quad \Sigma=\begin{bmatrix}0.2259 & 0.2259\end{bmatrix}$$ for which the log-likelihood was $L_2=-6.2123$.

However, $L_2 < L_1$! Even though I programmed the algorithm in python and got the same results, I know that $L_i$ ($i=0,1,2,\ldots)$ should be monotonically increasing, so, there has to be some mistake, but I can't catch it. Can anyone help?

Answer 1

Sorry, that this response is almost 2 years late. The OP has probably moved on, but in case anyone else encounters the same problem here is the solution.

First, note that I programmed the problem and then tried common errors until I was able to replicate your exact (erroneous) output:

initial parameters:
--------------------------------------------------
mu: (1, 3)
sigma: (1, 1)
pi: [0.5 0.5]

log_like: -4.389253938647963

iteration: 1
--------------------------------------------------
responsibilites: 
[[0.88079708 0.11920292]
 [0.5        0.5       ]
 [0.11920292 0.88079708]]

N_k: [1.5 1.5]
mu: [1.49227056 2.50772944]
sigma: [0.40887749 0.40887749]
pi: [0.5 0.5]

log_like: -3.675991903766752

iteration: 2
--------------------------------------------------
responsibilites: 
[[0.99770338 0.00229662]
 [0.5        0.5       ]
 [0.00229662 0.99770338]]

N_k: [1.5 1.5]
mu: [1.3363955 2.6636045]
sigma: [0.22629573 0.22629573]
pi: [0.5 0.5]

log_like: -6.19482298137464

The fix is simple. The $\Sigma$ returned by Bishop is the variance. But somewhere in your code you are using it as the standard deviation. Likely when you need to compute $N(\mathbf{x}_n| \mathbf{\mu}_k, \mathbf{\Sigma}_k)$ (as given by Bishop). Again, Bishop uses the notation that the normal distribution is parameterized by its mean and variance. But others like to think of the normal distribution as being specified by its mean and standard deviation. Both viewpoints are great, you just need to know which one you are using.

You mentioned that you coded this in Python. Are you using norm from scipy.stats? If so, this norm expects the standard deviation for the scale argument. So you need to take the square root of the $\Sigma$ returned by Bishop's equations before passing it as the scale argument to scipy.stats.norm.