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.
 A: 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!
A: 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.
