function [EE,AA]=Baum_Welch(observation,emmision,transmat,alpha1,beta1,EE,AA);
%%%%%%%%% Estimation    E%%%%%%%%%%
D=[5 8  11  14  17 20];
I=[2 4  7   10  13 16  19];
M=[1 3  6   9   12 15  18  21];
Q=6;
E=zeros(3*Q+3,20);
T=length(observation);
p=sum(alpha1(:,T));
   for k=1:Q+1  
           for a=1:20
              for i=1:T
                 if observation(i)==a 
                 if k==1
                     E(I(k),a)=E(I(k),a)+(1/p)*alpha1(I(k),i)*beta1(I(k),i);
                 else
                     E(M(k),a)=E(M(k),a)+(1/p)*alpha1(M(k),i)*beta1(M(k),i);
                     E(I(k),a)=E(I(k),a)+(1/p)*alpha1(I(k),i)*beta1(I(k),i);
                 end
                 end
              end
           end
   end
                  
 
    %%%%%%%%% Estimation   A%%%%%%%%%%
   A=zeros(3*Q+3,3*Q+3);
   emmision(M(Q+2),T+1)=1;
   for k=1:(Q+1)
            if k==1
            for i=1:T-1
             A(M(k),M(k+1))=A(M(k),M(k+1))+(1/p)*transmat(M(k),M(k+1))*emmision(M(k+1),observation(i))*beta1(M(k+1),i);
             A(I(k),M(k+1))=A(I(k),M(k+1))+(1/p)*alpha1(I(k),i)*transmat(I(k),M(k+1))*emmision(M(k+1),observation(i))*beta1(M(k+1),i);
             A(M(k),I(k))=A(M(k),I(k))+(1/p)*transmat(M(k),I(k))*emmision(I(k),observation(i))*beta1(I(k),i); 
             A(I(k),I(k))=A(I(k),I(k))+(1/p)*alpha1(I(k),i)*transmat(I(k),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1); 
             A(I(k),D(k))=A(I(k),D(k))+(1/p)*alpha1(I(k),i)*transmat(I(k),D(k))*beta1(D(k),i); 
            end
            observation=[21,observation];
            T=length(observation);
            A(M(k),D(k))=A(M(k),D(k))+(1/p)*transmat(M(k),D(k))*beta1(D(k),1);
            observation=observation(2:T);
            T=length(observation);
            elseif k==(Q+1)
              for i=1:T-1   
              A(M(k),M(k+1))=(1/p)*alpha1(M(k),T)*transmat(M(k),M(k+1));  
              A(M(k),I(k))=A(M(k),I(k))+(1/p)*alpha1(M(k),i)*transmat(M(k),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1);
              A(I(k),I(k))=A(I(k),I(k))+(1/p)*alpha1(I(k),i)*transmat(I(k),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1); 
              A(I(k),M(k+1))=(1/p)*alpha1(I(k),T)*transmat(I(k),M(k+1));
              A(D(k-1),M(k+1))=(1/p)*alpha1(D(k-1),T)*transmat(D(k-1),M(k+1));
             end
            else
             for i=1:T-1  
             A(M(k),M(k+1))=A(M(k),M(k+1))+(1/p)*alpha1(M(k),i)*transmat(M(k),M(k+1))*emmision(M(k+1),observation(i+1))*beta1(M(k+1),i+1);
             A(M(k),I(k))=A(M(k),I(k))+(1/p)*alpha1(M(k),i)*transmat(M(k),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1); 
             A(I(k),M(k+1))=A(I(k),M(k+1))+(1/p)*alpha1(I(k),i)*transmat(I(k),M(k+1))*emmision(M(k+1),observation(i+1))*beta1(M(k+1),i+1);
             A(I(k),I(k))=A(I(k),I(k))+(1/p)*alpha1(I(k),i)*transmat(I(k),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1); 
             A(D(k-1),M(k+1))=A(D(k-1),M(k+1))+(1/p)*alpha1(D(k-1),i)*transmat(D(k-1),M(k+1))*emmision(M(k+1),observation(i+1))*beta1(M(k+1),i+1);
             A(D(k-1),I(k))=A(D(k-1),I(k))+(1/p)*alpha1(D(k-1),i)*transmat(D(k-1),I(k))*emmision(I(k),observation(i+1))*beta1(I(k),i+1);
             end
             for i=1:T
             A(D(k-1),D(k))=A(D(k-1),D(k))+(1/p)*alpha1(D(k-1),i)*transmat(D(k-1),D(k))*beta1(D(k),i);
             A(M(k),D(k))=A(M(k),D(k))+(1/p)*alpha1(M(k),i)*transmat(M(k),D(k))*beta1(D(k),i);
             A(I(k),D(k))=A(I(k),D(k))+(1/p)*alpha1(I(k),i)*transmat(I(k),D(k))*beta1(D(k),i); 
           
             end
        end
   end
   A(M(Q+2),M(1))=1;
   AA=AA+A;
   EE=EE+E;
   
   
   
 
   