I am using the bootstrap algorithm to compute standard errors of the estimates of my normalmixEM output. I am not really sure if they are reliable?

My code is (data here):

# load package

B = 1000 # Number of bootstrap samples
mu1sample <- mu2sample <- sigma1sample <- sigma2sample <- lambdasample <- vector()

# Bootstrap
for(i in 1:B){
  subsample = sample(mydatatest,rep=T)
  normalmix <- normalmixEM(subsample, mu=c(-0.002294,0.002866),sigma=c(0.00836,0.02196), lambda=c(0.6746903,(1-0.6746903)),k=2, fast=FALSE, maxit=10000, epsilon = 1e-16, maxrestarts=1000)
  mu1sample[i]    = normalmix$mu[1]      # $
  mu2sample[i]    = normalmix$mu[2]      # $
  sigma1sample[i] = normalmix$sigma[1]   # $
  sigma2sample[i] = normalmix$sigma[2]   # $
  lambdasample[i] = normalmix$lambda[1]  # $
# standard errors


# show distribution of the bootstrap samples

This gives the following pictures:











EDIT: If you look at my variable, the mydatatest and use a KD to show the distribution with the following code


it looks like tt

2nd EDIT: I now included the mus and sigmas to be fixed. I updated the code and the pictures. Now again my question, what do you think about it?

  • 2
    $\begingroup$ What does "unreliable" mean here? Puzzling? I can't read your data but (for example) a bimodal bootstrap sample is very common when outliers are present: the outlier(s) will appear in some samples but not others. There are many relatives of this situation. Closer inspection of the raw data can often make much clearer what is going on. More to the point, it seems that you are fitting a mixture of normals. Is this exactly what you might expect with a mixture of normals? $\endgroup$ – Nick Cox Aug 10 '13 at 17:54
  • $\begingroup$ Sorry, you have to add an ".RData" after the filename to make it work. Yes I am fitting a mixture of normals, I want to have the standard errors of the estimates. I hope you can read my data now? $\endgroup$ – Jen Bohold Aug 10 '13 at 18:38
  • $\begingroup$ Too much like hard work for me if I have to read it into R first, not a language I use routinely. The real question is whether you take my point or can rebut it as not the right explanation. $\endgroup$ – Nick Cox Aug 10 '13 at 18:48
  • $\begingroup$ @NickCox Well, I am not 100% what you are aiming at? I am fitting a mixture of normals (two normals). So there is a market regime which represents the extreme "market scenario". I guess that's why the lambda estimate has two peaks? The results of the algorithm are not always the same if you run it without bootstrap. So sometimes the probability of pi (or lambda) (scenario 1) is e.g. 0.3 in the next time the vice versa way of 0.7. Sometimes even other numbers occur? I am not a professional on this topic, so maybe you can tell me what you need further to know? $\endgroup$ – Jen Bohold Aug 11 '13 at 8:58
  • $\begingroup$ @NickCox I added a kernel density of my data. Maybe not the typical data for a mixture of normals, but later on I do it with a different data set. And there it makes sense, because there is the big bell curve and thenn a small additional peak in the right tail, which I want to model with the second normal distribution. So the mixture will contain the big bell curve and the small additional peak. $\endgroup$ – Jen Bohold Aug 11 '13 at 9:03

As explained by Nick Cox and an anonymous user, what you think of as instability is just what the mixture models do: they don't care about labels unless you make it very clear that you know what your modes look like, roughly.

In terms of what you can do about fixing the labels where you need them to be, you would want to feed the full sample estimates of everything (both $\mu$s, both $\sigma$s, not just the $\lambda$ that you are feeding in now) as starting values. One can argue that this violates the spirit of maximum likelihood, but that may be the best you can do. If that does not really work, you may have to force even more information in, like insisting that $\mu_1 < \mu_2 - \delta$ and $\sigma_1 < \sigma_2 - \Delta$ and $\lambda > \frac12$. If normalmixEM() does not support that kind of cruelty to the parameter space, you would need to write your own likelihood with your own parameterization that accounts for such relations.

| cite | improve this answer | |
  • $\begingroup$ ok thanks for the answer and how can I fix the mus and sigmas and not only the lambda? What would the code look like? $\endgroup$ – Jen Bohold Aug 13 '13 at 8:28
  • $\begingroup$ I don't really know -- study the normalmixEM documentation. I wrote my own Gaussian mixture code in Stata some 15 years ago. Writing code to solve specific problems helps a lot in figuring things out about the model, in general. $\endgroup$ – StasK Aug 13 '13 at 15:14
  • $\begingroup$ I updated my answer according to your answer. What do you think about it now? $\endgroup$ – Jen Bohold Aug 14 '13 at 7:54
  • $\begingroup$ If sigma2 comes out bigger, swap the sigmas and the means you save from each run. $\endgroup$ – Nick Cox Aug 14 '13 at 8:00
  • $\begingroup$ @Jen, I think that's much better now. $\endgroup$ – StasK Aug 14 '13 at 12:47

My hunch is that your approach might not be reliable due to label switching, that is, each time you fit the mixture model, it's possible that the roles of the two normal distributions has been reversed.

That is, for different runs of the EM algorithm (mu1, sigma1) and (mu2, sigma2) might be switching roles.

It looks like the boot.se function provided by mixtools tries to account for this issue.

| cite | improve this answer | |
  • $\begingroup$ thanks for you answer, but the problem is, that boot.se leads to very different results every time I run the procedure? $\endgroup$ – Jen Bohold Aug 13 '13 at 8:29

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.