mixnorm <- function(y, mu0, mu1, var0, var1, prob, eps = 1/100000) {
    # y are assumed to be a mixture of two normal distributions
      new.params <- c(mu0, mu1, var0, var1, prob)
      err <- 1
      iter <- 1
      maxiter <- 4
      hist(y,probability=T,nclass=30)
      xvals <- seq(from=min(y),to=max(y),length=100)
      while(err > eps) {
       if (iter <= maxiter) {
         lines(xvals,prob*dnorm(xvals,mu1,sqrt(var1))+
             (1-prob)*dnorm(xvals,mu0,sqrt(var0)),lty=iter)
       }
      bayes <- (prob * dnorm(y, mu1, sqrt(var1)))/((prob * dnorm(y,
           mu1, sqrt(var1))) + ((1 - prob) * dnorm(y, mu0, sqrt(
           var0))))
       mu1 <- sum(bayes * y)/sum(bayes)
       mu0 <- sum((1 - bayes) * y)/sum((1 - bayes))
       var1 <- sum(bayes * (y - mu1)^2)/sum(bayes)
       var0 <- sum((1 - bayes) * (y - mu0)^2)/sum((1 - bayes))
       prob <- mean(bayes)
       old.params <- new.params
       new.params <- c(mu0, mu1, var0, var1, prob)
       err <- max(abs((old.params - new.params)/new.params))
       iter <- iter + 1
    }
  legend(4,.25,legend=c("iter 1","iter 2","iter 3","iter 4"),lty=1:4)
  return(list(mu = c(mu0, mu1), v = c(var0, var1), p = prob))
  }
    
  y <- c(rnorm(800,0,1),rnorm(200,3,3))
    
  vals <- mixnorm(y,-1,2,1,1.1,.7)