Expectation maximization and Gaussian mixtures - bad results

Question

I am supposed to find parameters of individual gaussians in a 1D mixture with a known number of components. I use my own implementation of EM algorithm; however, I am not able to find the right solution.

Data

Number of points: 5057
Number of components: 5

Problem

Everytime I launch the algorithm, all resulting component's means are values from 161 to 163 (That is between the first two and largest peaks).

Here are my results from individual steps of one run of EM: output
Here is my implementation of EM in C++: EM

Update

Pat, you were right, thank you for your help.

1. I forgot to multiply each item of product by posterior probability
2. I had mistake in log likehood function

Answer 1

There's a bug in your code, is the simple answer. And honestly, I'm not going to spend too long slogging through someone else's custom c++ implementation to try and figure out where (even if yours is laid out in a commendably clear fashion).

Assuming it's not just a typo somewhere, things you might want to look at:

• Have you checked the log likelihood is increasing with every iteration? Right now you have a criterion function which looks similar to the log likelihood, but not quite the same, and it is decreasing each iteration. Unless there's a minus sign somewhere in there that I've missed, that's a bad sign. Always always always check this with EM algorithms.
• QVector<double> v = MathUtils::sub(x, g->getMean()); g->setSigma(sqrt(MathUtils::dot(v, v) / sum)); Are you sure this is giving the correct value for sigma? In particular, I see no point where you weight the vector v by the posterior probabilities in probs (Compare this to the previous line, which uses MathUtils::dot(x, probs) / sum when calculating the mean). If I had to bet, I'd say this is where the bug may be - I'd guess it should be something like g->setSigma(sqrt(MathUtils::dot(v, MathUtils::dot(v,probs)) / sum));. EDIT: Actually, now I look at it, my guess at what the code should look like is wrong - that dot product will return a scalar, and you don't want that, you want a vector. However, by the sounds of things, you've fixed the issue yourself
• Are you certain hs.set(i, j, bayes->eval(x.at(i))); is giving you the posterior probabilities? I can't figure out what this bayes object is exactly, but make sure that it's giving values in the range 0 to 1, and for each sample they sum to 1.

Other than that, I've not got any suggestions. Happy bug hunting!