# The first example in Dempster et al.

y <- c(125,18,20,34)
log.lik <- function(u) y[1]*log(2+u)+(y[2]+y[3])*log(1-u)+y[4]*log(u)
u <- seq(0.001,0.999,length=100)
plot(u,log.lik(u),type="l")

?optimize
optimize(log.lik,c(0,1),maximum=TRUE)
u.opt <- optimize(log.lik,c(0,1),maximum=TRUE)$maximum

u.em <- 0.5
for(i in 1:8)
{
    x <- y[1]*u.em/(2+u.em)
    u.em <- (x+y[4])/(x+sum(y[2:4]))
}
u.em

u.em <- 0.5
U.em <- rep(u.em,8)
for(i in 1:8)
{
    x <- y[1]*u.em/(2+u.em)
    u.em <- (x+y[4])/(x+sum(y[2:4]))
    U.em[i] <- u.em
}
cbind(U.em,U.em-u.opt,(U.em-u.opt)/(c(0.5,U.em[1:7])-u.opt))




# An Example of Normal Mixture
N <- 10
alpha <- c(0.5,0.25,0.25)
mu <- (1:3)*5
mu.y <- sample(mu,N,replace=TRUE,prob=alpha)
y <- rnorm(N,mu.y)
mu.y
y

alpha0 <- c(0.1,0.4,0.5)
mu0 <- 1:3


k<-length(alpha0)
n<-length(y)
L<-exp(-(outer(y,mu0,"-"))^2/2)
P<-sweep(L,2,alpha0,"*")
sP<-rowSums(P)
P<-sweep(P,1,sP,"/")
alpha<-apply(P,2,mean)
mu<-colMeans(sweep(P,1,y,"*"))/alpha

tol<-max(abs(c(mu-mu0,alpha-alpha0)))
mu0<-mu
alpha0<-alpha

EMmix<-function(alpha0,mu0,y,Tol=10^(-5),MaxIter=10^3)
{
    iter<-1
    tol<-1
    while((iter<=MaxIter) & (tol>Tol))
    {
        L<-exp(-(outer(y,mu0,"-"))^2/2)
        P<-sweep(L,2,alpha0,"*")
        sP<-rowSums(P)
        P<-sweep(P,1,sP,"/")
        alpha<-colMeans(P)
        mu<-colMeans(sweep(P,1,y,"*"))/alpha
        tol<-max(abs(c(mu-mu0,alpha-alpha0)))
        mu0<-mu
        alpha0<-alpha
        iter<-iter+1
    }
    if(iter > MaxIter) warning("Maximal number of iterations exceeded")
    return(list(alpha=alpha,mu=mu,iter=iter-1))
}
alpha0 <- c(0.1,0.4,0.5)
mu0 <- 1:3
EMmix(alpha0,mu0,y)

N <- 1000
alpha <- c(0.5,0.25,0.25)
mu <- (1:3)*5
mu.y <- sample(mu,N,replace=TRUE,prob=alpha)
y <- rnorm(N,mu.y)
alpha0 <- c(0.1,0.4,0.5)
mu0 <- c(5,7,2)
EMmix(alpha0,mu0,y)

alpha0 <- rep(1,3)/3
mu0 <- rep(mean(y),3)
EMmix(alpha0,mu0,y)