EM.Anae=function() {
data.work=options()$data.work
m=dim(data.work)[[2]]-2
no.iter=options()$no.iter
K=options()$K
pi1=options()$pi1
credmatrix=options()$credmatrix

## the last column is fully determined by the first three columns

continue=T
value1=-1
count1=0
l1=1

while (continue) {
count2=0
continue1=T
while (continue1) {
    count2=count2+1
    print(paste("iterative updating likelihood function ", count2, sep=""))
    for (j in 1:m) {
        options(R=j)
        options(credmatrix=credmatrix)
        cred.temp=credmatrix[((j-1)*K*(K-1)+1):((j-1)*K*(K-1)+K*(K-1))]
        parameter=optim(cred.temp,likelihood.Anae, method = c("Nelder-Mead"), list(maxit=600, abstol=1.e-6,reltol=1.e-6, alpha=1, beta=0.5, gamma=0.5))
        credmatrix[((j-1)*K*(K-1)+1):((j-1)*K*(K-1)+K*(K-1))]=parameter$par
        print(parameter$value)
        if (j ==1) l=parameter$value;
    }
    
    parameter=optim(credmatrix,likelihood.O.Anae, method = c("Nelder-Mead"), list(maxit=100, abstol=1.e-6,reltol=1.e-6, alpha=1, beta=0.5, gamma=0.5))
    credmatrix=parameter$par
    l1=parameter$value;
    continue1=F;
    if( abs(l-l1)>1) {continue1=T;}
}

#credmatrix1=credmatrix;options(credmatrix=credmatrix1)
#pi2=pi1;options(pi1=pi1)

if(1) {
## update conditional probability matrix
    n1=dim(data.work)[[1]]
    m1=dim(data.work)[[2]]
## p(k|x1)=p(x1|k)p(k)/sum(p(x1|k)p(k))      # p(k) are the prior prob 1/K
    for (ii in 1:n1) {
            temp=rep(0,K)
            for(jj in 1:m1) {
                kk=ifelse(jj>3,jj-2,1 )
                cred.temp=credmatrix[((kk-1)*K*(K-1)+1):((kk-1)*K*(K-1)+K*(K-1))];
                cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
                cred.temp=cbind(cred.temp, 1-rowSums(cred.temp))
                obs=data.work[ii,jj]
                a=cred.temp[, obs] *pi1[ii,]
                temp=a/sum(a)+temp
            }
##           temp=temp*pi1[ii, ]
### Expectation
        pi1[ii,]=temp/sum(temp)   
    }
}
    value2=parameter$value
    print(paste(value2, ", diff= ", round(abs(value1-value2), 2)))
    if( abs(value1-value2)<0.001  ) {continue=F; print(paste("Iteratinon Converges at ", count1))}
    else {value1=value2}
    #if (count1==no.iter) {print("maximum iteratoin is reached and program is not converged "); continue=F}
    options(value=value2)
    options(pi1=pi1)
} 

}