# Algorithms for 1D Gaussian mixture with equal variance and noise cluster

I would like to fit a Gaussian mixture model to some data. The data is 1D and I want to constrain all the Gaussians to have equal variance. I would also like to have a uniform background noise cluster to pick up points that are not in a real cluster.

I don't know how many clusters there should be, and this can vary over a large range for different data sets (tens to thousands are possible), so ideally this should be selected automatically. It might also be useful to put some lower threshold on the mixing coefficient for each cluster to avoid fitting very tiny ones to the background noise.

Depending on the data set there will be from hundreds to many thousands of data points.

What would be a good way to pursue this problem? I am a bit overwhelmed by the possible options (EM, variational bayes, infinite GMM, dirichlet processes etc.), or if it is even possible.

As I understand, there are three questions in one, I will try to suggest some answers for all of them.

1. How to estimate the number of clusters?

In 1D I ones tried to estimate the number of clusters by density estimation. First you try to estimate density (for example, see ksmooth in R) with automatically established parameters. Then you can approximate the number of clusters by just evaluating the number of local maximums of the density. That would be a good number to start with, after that you can do cross-valudation (for the final method) to get the more exact number.

2. How to fit the models with equal variance?

That can be easily incorporated into EM-procedure. During M-step you estimate the mean $\mu_k$ as usual, but the variance $\sigma_k$ like this: $\sigma_k = \sigma = \frac{1}{N}\sum_{i=1}^N (x_i - \mu_{k(x_i)})$, where $k(x_i)$ return an estimated class of the object $x_i$.

3. How to deal with uniform background noise?

If it's really uniformly distributed, it should not greatly affect the EM-procedure. And after all, you can get rid of it, just marking $x_i$ belonging to noisy class if the probability of this object belonging to any class $k$ is less then some threshold. You can think of incorporating this step inside EM-algorithm, but I am not sure how it would work.

Hope some of the suggestions above would help!

Dmitry's answer is mostly correct, but there's a major problem with the second point: estimating variance under the assumption they're all equal. Firstly, it should be $\sigma^2$, not $\sigma$. Secondly, we need to find the L2 distance, not L1, so square the difference term: $\sigma^2=\frac{1}{N}\sum_{i=1}^N(x_i-\mu_{k})^2$. Thirdly, this does not even actually work. This unfortunately does not give the optimal variance. I ran a bunch of simulations, and it does not converge to the local minimum. The problem is that if you try to solve for the variance by taking the derivative of the PDF and setting it to 0, you don't get a closed form for $\sigma^2$. Here's the PDF:$$L(\mu,\omega,\sigma^2|x)=\prod_{i=1}^N\left(\sum_{k=1}^K\omega_k\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x_i-\mu_k)^2}{2\sigma^2}}\right),$$ so the log likelihood becomes:$$l(\mu,\omega,\sigma^2|x)=\sum_{i=1}^Nlog\left(\sum_{k=1}^K\omega_k\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x_i-\mu_k)^2}{2\sigma^2}}\right).$$ If we take the derivative and set it to 0, we get: $$\sigma^2=\frac{\sum_{i=1}^N\frac{\sum_{k=1}^K2\omega_k(x_i-\mu_k)^2N(\mu_k,\sigma)}{\sum_{k=1}^K\omega_kN(\mu_k,\sigma)}}{N},$$ which as you can see is not a closed form. Very unfortunate. There's not much you can do besides trying to solve numerically. There are many methods to do this, but depending on your problem you may want to do one over another.

My only hope is that I'm wrong, but I don't think I am. I went through the derivation several times, and found no errors. If you have already solved this problem, please let me know!

VERY IMPORTANT EDIT:

After some more research, I found our answer here: http://lasa.epfl.ch/teaching/lectures/ML_Phd/Notes/GP-GMM.pdf

I'll paraphrase what the paper says. While the above formulation is correct, there is another way to describe the model. If we consider each point $x_i$ to be generated by gaussian $k$, we can think of the total likelihood to be the product of the likelihoods of the $K$ gaussians, as they are independent. Of course, when we do the E step of the EM algorithm, we come up with a posterior probability for each point, $P(y=k|x_i)$, where $y=1,2,...,K$ is the gaussian that generated point $x_i$. So the likelihood for any gaussian will be: $$L(\mu_k,\sigma^2|x,y=k)=\prod_{i=1}^N P(x_i)\cdot P(y=k|x_i), \\ L(\mu_k,\sigma^2|x,y=k)=\prod_{i=1}^N\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x_i-\mu_k)^2}{2\sigma^2}}\cdot P(y=k|x_i),$$ SO: the full likelihood (assuming the three Gaussians are independent) is just the product:$$L(\mu_k,\sigma^2|x)=\prod_{k=1}^K\prod_{i=1}^N\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x_i-\mu_k)^2}{2\sigma^2}}\cdot P(y=k|x_i) ,\\log L(\mu_k,\sigma^2|x) =l(\mu_k,\sigma^2|x) =\sum_{k=1}^K\sum_{i=1}^Nlog\left(\frac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{-(x_i-\mu_k)^2}{2\sigma^2}}\right)\cdot P(y=k|x_i).$$ This makes it super easy to optimize:$$\mu_k=\frac{\sum_{i=1}^Nx_ip(k|x_i)}{N},\\ \omega_k=\frac{\sum_{i=1}^Np(k|x_i)}{N},$$ and the variance (total variance across all gaussians): $$\sigma^2=\frac{\sum_{k=1}^K\sum_{i=1}^N (x_i-\mu_k)^2\cdot P(y=k|x_i)}{N}.$$ So that's it! I hope it comes in useful for someone.

• I never got very far with the problem - I tried tweaking some implementations I found online but I can't remember the details now. But the likelihood you have there is completely general (ie for any cluster assignment). In the EM procedure you have a canditate point-cluster assignment (given by k(xi) in Dmitry's answer) so I think it should be L(mu, omega, sigma| x, k) and then the second sum over k is not necessary (just uses k(xi) as a function of x) and then I think you would get the normal ML estimate for sigma? – thrope Sep 22 '14 at 7:35
• Are you saying the second sum is unnecessary because you can classify each point to be in a particular cluster? If so, this runs into issues, as you're essentially doing a K-means then, instead of adding a probabilistic value to each point belonging in a given cluster... – ArturJ Oct 7 '14 at 18:43
• That's what my lab mate says -- to essentially maximize over all points, conditioning on each cluster k. What I think ends up happening then, is you end up taking a sort of weighted average of the three clusters' variances, which unfortunately does not give a maximal likelihood. When you take the partial of $l(\mu,\omega,\sigma^2|x,k)$ with respect to $\sigma^2$, you still have that conditioning on the cluster k, while you should have the general solution. Anyway, I tried conditioning on k, and the likelihood doesn't decrease monotonically, suggesting an error. – ArturJ Oct 7 '14 at 19:00
• I don't really understand the point about k-means. I do think the second sum is unnecessary, because I think in the M step you are calculating the mean and standard deviation for the specific candidate cluster assignment (so xi and k(xi)). These will then be used in the E-step to determine the next candidate assignments. So I think the log-likelihood should be the same as in your second expression, but with the sum over k removed and k's replaced with k(xi). – thrope Oct 7 '14 at 20:58