Hi all,

I'm trying to create an em algorithm for mixture problems but my code has a mistake that I can't find.

Y<-c(5,9,8,4,7)

ln<-function(theta,Y){

-sum(log(theta[1]*dbinom(Y,10,theta[2])+(1-theta[1])*dbinom(Y,10,theta[3])))

}

#EM algorithm

emstep<-function(y,theta){

EZ<-theta[1]*dbinom(Y,10,theta[2])/(theta[1]*dbinom(Y,10,theta[2])+(1-theta[1])*dbinom(Y,10,theta[3]))
theta[1]<-mean(EZ)
theta[2]<-sum(EZ*Y)/sum(EZ)

theta[3]<-sum((1-EZ)*Y)/sum(1-EZ)

theta

}

emiteration<-function(Y,theta,k){

for(i in 1:k){

paliotheta<-theta

theta<-emstep(Y,theta)

if((paliotheta[1]-theta[1])^2+(paliotheta[2]-theta[2])^2+(paliotheta[3]-theta[3])^2<=1e-10) {

break

}

}

theta

}

#starting values

theta<-c(0.5,0.4,0.3)

#1000 iterations of EM algorithm

theta<-emiteration(Y,theta,1000)

theta

#check for convergence

theta<-emstep(Y,theta)

theta

hist(Y)

I'd really appreciate any help.

Thanks!