EM algorithm with a component for outliers i have a vector of measurements from one to three classes, which can be modeled by gaussian distributions. There are some outliers in the data. I use the EM algorithm to learn the parameters of the components. One property of the algorithm is, that the incomplete likelihood of the data increases monotonically. My code works fine, but the parameters are heavily biased due to the outliers. Now i want an extra component, a uniform distribution, to model the outliers. And now the monotonicity is no longer guaranteed. As far as i know the new component only changes the E step of the algorithm. Does anyone know what my mistake is?
library(gtools)

#create 200 data points from two components plus one outlier
y = c(rnorm(n=100,mean=-4,sd=.5),rnorm(n=100,mean=0,sd=.5),4)
stripchart(y,method='jitter',pch=4)


em = function(x,maxit){

K = 2 #number of components
n = length(x) #size of dataset

#start values
p = 1 / (max(x)-min(x)) #probability for the uniform distribution
pi = rdirichlet(n=1,alpha=c(rep(5,times=K),1))[1,] # mixing probabilities
mu = sample(x,size=K) #vector of means
sigma = .7 #standard deviation
gamma = 0 # 

#maximization
m_step = function(){
    pi = apply(gamma,2,sum)
    mu = gamma[,1:K] * x
    mu = apply(mu,2,sum)
    mu = mu / pi[1:K]
    pi = pi / n
    sigma = t(sapply(x,`-`,mu))
    sigma = gamma[,1:K] * sigma^2
    sigma = sqrt(sum(sigma) / n)
    pi <<- pi; mu <<- mu; sigma <<- sigma
    return()
}

#expectation
e_step = function(){
    gamma = sapply(1:K,function(k){ dnorm(x,mean=mu[k],sd=sigma)},USE.NAMES=FALSE)
    gamma = cbind(gamma,p)
    gamma = t(t(gamma)*pi)
    tmp = apply(gamma,1,sum)
    gamma <<- gamma / tmp
    loglik = sum(log(tmp))
}

loglik = rep(NA,times=maxit)

loglik[1] = e_step()
m_step()

j = 1
gain = Inf
while(j <= maxit && gain > 1e-8){
    loglik[j+1] = e_step()
    gain = loglik[j+1] - loglik[j]
    m_step()
    j = j+1
}

if(any(diff(loglik[1:j])<0)){
    print('Error: log likelihood not monotonic')
    print(loglik[1:j])
}

return(list(pi=pi,mu=mu,sigma=sigma,loglik=loglik[j]))
}

print(em(y,200))

I found a solution. I introduced a new hidden binary variable with the meaning outlier/no outlier and now it works. I will try to write down the math and post it later.
 A: I think I may have an explaining of this phenomenon, although this would require some experiments to be sure this is the true reason... anyway I hope this helps.
After introducing this "outlier component", your likelihood landscape becomes very messy : you have plenty of spikes (on the border of the parameter space) going to an infinite likelihood, when you take e.g. $\mu_1 = {}$ one of the data points, says $x_1$, and you make $\sigma$ go to 0 (and $\mu_2$ is wherever you want). 
In that case, you will have one data point ($x_1$) contributing to the likelihood with a very huge value (as $\sigma$ is very small), and all other data points will be classified as outliers (except possibly one point at $\mu_2$...) and will contribute to the likelihood by a value $\simeq p$ (using the notations in your script). It will happen that from time to time (depending on your starting values) the EM is going to converge such a point. 
I suspect (this is what I haven't tested) that when this happens, some of the computations become unreliable, because of rounding problems in the steps where you work with likelihoods or densities and not with log-likelihoods. 
You may model your outliers by a normal component with standard deviation, say, $10 \sigma$, I think this will (partly) avoid this problem...
