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(EZY)/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!