eq<-c(13, 14, 8, 10, 16, 26, 32, 27, 18, 32, 36, 24, 22, 23, 22, 18, 25, 21, 21, 14,
8, 11, 14, 23, 18, 17, 19, 20, 22, 19, 13, 26, 13, 14, 22, 24, 21, 22, 26, 21,
23, 24, 27, 41, 31, 27, 35, 26, 28, 36, 39, 21, 17, 22, 17, 19, 15, 34, 10, 15,
22, 18, 15, 20, 15, 22, 19, 16, 30, 27, 29, 23, 20, 16, 21, 21, 25, 16, 18, 15,
18, 14, 10, 15, 8, 15, 6,11, 8, 7, 18, 16, 13, 12, 13, 20, 15, 16, 12, 18,
15, 16, 13, 15, 16, 11, 11)


l0<--1000000
EMmix<-function(data,p,lambda1,lambda2){
n<-length(data)
z1<-rep(0,n)
test<-3
while (test!=10) {
for (i in 1:n){
z1[i]<-p*dpois(data[i],lambda1)/(p*dpois(data[i],lambda1)+(1-p)*dpois(data[i],lambda2))
}
z2<-1-z1
p.new<-sum(z1)/n
lambda1.new<-sum(data*z1)/sum(z1)
lambda2.new<-sum(data*z2)/sum(z2)
lnew<-sum(z1)*log(p.new)+sum(z1*dpois(data,lambda1.new,log=TRUE))+sum(z2)*log(1-p.new)+sum(z2*dpois(data,lambda2.new,log=TRUE))
if (abs((lnew-l0)/lnew)<=10^(-6)) {
                 test<-10} else {
                     print(c(lnew-l0,lambda1.new,lambda2.new))
                     lambda1<-lambda1.new
                     lambda2<-lambda2.new
                     p<-p.new
                     l0<-lnew
                      }
}  
return(list(lambda1.new=lambda1.new,lambda2.new=lambda2.new,p.new=p.new,z1=z1))
}