3
$\begingroup$

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.

$\endgroup$

1 Answer 1

2
$\begingroup$

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...

$\endgroup$
2
  • $\begingroup$ It is quite often the problem, that the outlier component has a large share of the data. But even when this is not the case, the problem occurs. And the difference is to big to be numerical only. Here for example the log likelihood from one run -583.3709 -420.7938 -407.7895 -416.2775 $\endgroup$
    – kalbfred
    Nov 14, 2013 at 13:42
  • $\begingroup$ This is very surprising. There might be something in the proof of monotonicity of the EM which fails in this situation (the uniform component is very special, as it has no parameter to be estimated in the M step) but I don’t see what. I don’t have enough time now to dig this, but if no one answers I promise to give it a more serious try. $\endgroup$
    – Elvis
    Nov 14, 2013 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.