0
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Momo Mar 30 '13 at 11:25
  • $\begingroup$ yes, I just compiled/run the SAME code several times, this gives different outputs! So I did NOT change any values @Momo ! $\endgroup$ – Stat Tistician Mar 30 '13 at 12:23
  • $\begingroup$ 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. $\endgroup$ – Momo Mar 30 '13 at 12:58
  • $\begingroup$ @Momo I did what you said and what should I do now? $\endgroup$ – Stat Tistician Mar 30 '13 at 13:45
  • $\begingroup$ 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. $\endgroup$ – Momo Mar 30 '13 at 15:37
2
$\begingroup$

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")
hist(dat,probability=T,add=T)
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ this is a real improvement and I will accept your answer, but one last thing: Sometimes I still get values for pi like [1] 0.9818572 so it is not really perfect, should I care about this or ignore this? $\endgroup$ – Stat Tistician Mar 30 '13 at 14:33
  • 1
    $\begingroup$ Are these results obtained with maixt=10000? I get always the same result if I specify this number of iterations. $\endgroup$ – Triple roast Mar 30 '13 at 14:37
  • 1
    $\begingroup$ 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. $\endgroup$ – Triple roast Mar 30 '13 at 14:48
  • 1
    $\begingroup$ Good. The only issue with the generalized hyperbolic distribution is that it has different tail behaviour in each direction while other skew-t distributions have the same tail behaviour. $\endgroup$ – Triple roast Mar 30 '13 at 14:51
  • 1
    $\begingroup$ Yeah, that is the problem with Guassian mixtures. The only way to recognise the best choice is by conducting a model comparison via some technique. Numerical methods are likely to get stuck in local maxima. Mixtures of distributions, although flexible, are more appropriate when the features of the data come from a mixing of several populations in the same data set. $\endgroup$ – Triple roast Mar 30 '13 at 14:54

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.