Stat 514: Spring 2010 

                                 Homework  7:



Use the below R code to report the final estimates of lamda, mu1 and mu2. Also give the total number of iterations to convergence.


#########################################################################
# below is a function to solve a normal mixture problem by EM algorithm #
#########################################################################

emnormal <- function(x,lambda,mu1,mu2,sigma=15,epsi=10^(-8),maxit=1000){
    diff <- 1
    iter <- 0
    oldlog <- sum(log(lambda*dnorm(x,mu1,sigma)+(1-lambda)*dnorm(x,mu2,sigma)))
    while(diff>epsi & iter<maxit){
        z <- lambda*dnorm(x,mu1,sigma)/(lambda*dnorm(x,mu1,sigma)+
        (1-lambda)*dnorm(x,mu2,sigma))
        lambda.new <- mean(z)
        mu1.new <- sum(z*x)/sum(z)
        mu2.new <- sum((1-z)*x)/sum(1-z)
        sigma.new=15
        #sigma.new <- sqrt((sum(z*(x-mu1.new)^2)+sum((1-z)*(x-mu2.new)^2))/length(x))
        lambda <- lambda.new
        mu1 <- mu1.new
        mu2 <- mu2.new
        sigma <- sigma.new
        newlog <- sum(log(lambda*dnorm(x,mu1,sigma)+(1-lambda)*dnorm(x,mu2,sigma)))
        diff <- newlog-oldlog
        oldlog <- newlog
        iter <- iter+1
    }
    if(iter==maxit){cat("warning! Not convergence","\n")}
    cat("number of iterations=",iter,"\n")
    list(posteri=z,lambda=lambda,mu1=mu1,mu2=mu2,sigma=sigma,logvalue=newlog)
}

x=c(105,98,114,104,86,120,100,124,87,99,98,
    138,99,112,114,76,114,111,100,126,100,95,
    102,105,80,115,81,93,126,87,104,119,94,
    139,96,105,84,77,114,122,120,106,104,143,
    120,93,102,84,136,71)
    
emnormal(x,.5,110,90)

#########################################################################
#   the end of the program                                              #
#########################################################################
# explanations                                                          #
# x: the data set
# lambda: the initial value of the mixing proportion, typically we choose .5
# mu1,mu2, the initial values for the two component's means
# sigma: the initial value for the variance NOTE:not standard deviation
# epsi, maxit: the controls
#########################################################################