# r: normalmixEM not consistent?

I have the following data: http://s000.tinyupload.com/?file_id=00083355432555420222

I want to fit a mixture density of two normal distributions.

I have the formula: \begin{align} f(l)=\pi \phi(l;\mu_1,\sigma^2_1)+(1-\pi)\phi(l;\mu_2,\sigma^2_2) \end{align}

my R code is:

normalmix<-normalmixEM(dat,k=2,fast=TRUE)

pi<-normalmix$lambda[1] mu1<-normalmix$mu[1]
mu2<-normalmix$mu[2] sigma1<-normalmix$sigma[1]
sigma2<-normalmix$sigma[2]  Now I have the problem, that the output is not consistent, i.e. every time I run the code, I get different outputs! And they are very different, no small differences, which could be due to the precision of the numerical procedures. E.g. sometimes for pi I get [1] 0.2653939  or [1] 0.3318069  I already recognized, that sometimes the numbering is changed, so the pi of 0.7 would be equal to a pi of 0.3. Okay, I got this, I don't know why the R procedures does this, but this would not be a problem. But the problem is, that the outputs are way to different, sometimes I even get an error message (german): Fehler in while (dl > eps && iter < maxit) { : Fehlender Wert, wo TRUE/FALSE nötig ist Also, the number of iterations is very different, from 29 up to 1000 ...... edit: I now tried to set the initial values manually.... In the Rmanual it says: Lambda: "Initial value of mixing proportions. Automatically repeated as necessary to pro- duce a vector of length k , then normalized to sum to 1" I started to try it with giving different values of lambda: 1. normalmix<-normalmixEM(dat,lambda=c(0.5,0.5),k=2,fast=TRUE) normalmix$loglik
normalmix$mu normalmix$sigma


which gives a loglik of 617.2996 and the parameters:

> normalmix$mu [1] 0.0003769442 -0.0008892282 > normalmix$sigma
[1] 0.02814997 0.01232638


2.

normalmix<-normalmixEM(dat,lambda=c(0.8,0.2),k=2,fast=TRUE)
normalmix$loglik normalmix$mu
normalmix$sigma  which also gives a loglik of 617.2996 ? and the parameters > normalmix$mu
[1]  0.0003767148 -0.0008892515
> normalmix$sigma [1] 0.02814819 0.01232488  mh, I don' know.... then, I tried it with different values of$\mu$: 1. normalmix<-normalmixEM(dat,mu=c(-0.001,0.002),k=2,fast=TRUE) normalmix$loglik
normalmix$mu normalmix$sigma


which also leads to a loglik of 617.2996 and the parameter output:

> normalmix$mu [1] 0.0003767138 -0.0008892516 > normalmix$sigma
[1] 0.02814818 0.01232487


2.

normalmix<-normalmixEM(dat,mu=c(0.02,0.01),k=2,fast=TRUE)
normalmix$loglik normalmix$mu
normalmix$sigma  which again! leads to a loglik of 617.2996 and the parameter output: > normalmix$mu
[1]  0.0003769462 -0.0008892280
> normalmix$sigma [1] 0.02814999 0.01232639 >  what should I do now? • This could habe to do with local minima for the mixture likelihood. Have you tried running the code with the same starting values every time? Does the behaviour still occur then? Imo error message is too unspecific to to be diagnosed from here. You can use traceback() to see where the error occurred. Also, possible bugs might better be reported to the package authors. – Momo Mar 30 '13 at 11:25 • yes, I just compiled/run the SAME code several times, this gives different outputs! So I did NOT change any values @Momo ! – Stat Tistician Mar 30 '13 at 12:23 • If you have not passed ANY values as starting values, the values will be RANDOMly generated. From the docs: "If NULL, then the initial value is randomly generated from a normal distribution determined by binning the data." So try passing the SAME starting values for mu and sigma a couple of times and see if it is still "inconsistent". If not, the starting values are probably the reason and the algorithm is stuck in local maxima. If so, run the algorithm with different starting values and take the solution that had the highest likelihood over all runs. – Momo Mar 30 '13 at 12:58 • @Momo I did what you said and what should I do now? – Stat Tistician Mar 30 '13 at 13:45 • Well, you now showed that the results are not inconsistent... So its just that the algorithm did not numerically find the optimum before termination or local maxima. – Momo Mar 30 '13 at 15:37 ## 1 Answer The problem is that the default number of iterations (maxit=1000) is not enough to guarantee the convergence of the EM algorithm in general and you have to fiddle a bit with this (and sometimes other) feature. It seems like maixt=10000 produces stable results. See the following code (run it several times to convince yourself). normalmix<-normalmixEM(dat,k=2,fast=TRUE,maxit=10000) pi<-normalmix$lambda[1]
mu1<-normalmix$mu[1] mu2<-normalmix$mu[2]
sigma1<-normalmix$sigma[1] sigma2<-normalmix$sigma[2]

c(pi,mu1,mu2,sigma1,sigma2)


Now, check the fit of the estimated mixture.

mix2 <- function(x) return( pi*dnorm(x,mu1,sigma1)+(1-pi)*dnorm(x,mu2,sigma2))

vec <- seq(-0.05,0.05,0.0001)

plot(vec,mix2(vec),type="l",col="red")

• Are these results obtained with maixt=10000? I get always the same result if I specify this number of iterations. – Triple roast Mar 30 '13 at 14:37
• I think you have to increase the precision, the number of initial points and the number of iterations. The rest is just related to the non-uniqueness of the estimation of MLE in a Gaussian mixture (as you mention). Check normalmix<-normalmixEM(dat,k=2,fast=FALSE,maxit=10000,epsilon = 1e-16,maxrestarts=1000). I would rather model these data using a skew-t in order to avoid local maxima. – Triple roast Mar 30 '13 at 14:48