function [beta1]=compute_beta(observation,emmision,transmat)
D=[5 8  11  14  17 20];
I=[2 4  7   10  13 16  19];
M=[1 3  6   9   12 15  18  21];
T=length(observation);
Q=6;
beta1=zeros((3*Q)+3,T);
beta1(M(Q+2),T+1)=1;
beta1(M(Q+1),T)=transmat(M(Q+1),M(Q+2));
beta1(I(Q+1),T)=transmat(I(Q+1),M(Q+2));
beta1(D(Q),T)=transmat(D(Q),M(Q+2));
for  k=Q:-1:1
    if k==1
       beta1(I(k),T)=beta1(D(k),T)*transmat(I(k),D(k));
    else
     beta1(D(k-1),T)=beta1(D(k),T)*transmat(D(k-1),D(k));  
     beta1(M(k),T)=beta1(D(k),T)*transmat(M(k),D(k));
     beta1(I(k),T)=beta1(D(k),T)*transmat(I(k),D(k));
    end
end

for i=T-1:-1:1
     for k=(Q+1):-1:1
          if k==Q+1
               beta1(I(k),i)= beta1(I(k),i+1)*transmat(I(k),I(k))*emmision(I(k),observation(i+1));
               beta1(M(k),i)=beta1(I(k),i+1)*transmat(M(k),I(k))*emmision(I(k),observation(i+1));  
               beta1(D(k-1),i)= beta1(I(k),i+1)*transmat(D(k-1),I(k))*emmision(I(k),observation(i+1));
          elseif k==1
              beta1(I(k),i)= beta1(M(k+1),i+1)*transmat(I(k),M(k+1))*emmision(M(k+1),observation(i+1))+ beta1(I(k),i+1)*transmat(I(k),I(k))*emmision(I(k),observation(i+1))+beta1(D(k),i)*transmat(I(k),D(k));
          else
               beta1(D(k-1),i)=beta1(M(k+1),i+1)*transmat(D(k-1),M(k+1))*emmision(M(k+1),observation(i+1))+ beta1(I(k),i+1)*transmat(D(k-1),I(k))*emmision(I(k),observation(i+1))+beta1(D(k),i)*transmat(D(k-1),D(k));
               beta1(I(k),i)=beta1(M(k+1),i+1)*transmat(I(k),M(k+1))*emmision(M(k+1),observation(i+1))+ beta1(I(k),i+1)*transmat(I(k),I(k))*emmision(I(k),observation(i+1))+beta1(D(k),i)*transmat(I(k),D(k));
               beta1(M(k),i)=beta1(M(k+1),i+1)*transmat(M(k),M(k+1))*emmision(M(k+1),observation(i+1))+ beta1(I(k),i+1)*transmat(M(k),I(k))*emmision(I(k),observation(i+1))+beta1(D(k),i)*transmat(M(k),D(k));
          end
     end
end
beta1=beta1(:,1:T);
