emmision2=emmision1;
emmision3=emmision2;
emmision4=emmision3;
transmat2=transmat1;
transmat3=transmat2;
transmat4=transmat3;
scoreml2=zeros(1,1000);
for i=1001:2000
index1=[];
index2=[];
index3=[];
index4=[];
p1=[];
p2=[];
p3=[];
p4=[];
k1=1;
k2=1;
k3=1;
k4=1;
  observation=all_ungap_sequences(i,find(all_ungap_sequences(i,:)>0));
  T=length(observation);
  [alpha1]=compute_alpha(observation,emmision1,transmat1);
  p1(1)=sum(alpha1(:,T));
  index1=[1];
  [alpha1]=compute_alpha(observation,emmision2,transmat2);
  p2(1)=sum(alpha1(:,T));
  index2=[1];
  [alpha1]=compute_alpha(observation,emmision3,transmat3);
  p3(1)=sum(alpha1(:,T));
  index3=[1];
  [alpha1]=compute_alpha(observation,emmision4,transmat4);
  p4(1)=sum(alpha1(:,T));
  index4=[1];
 for k=2:2000  
  p(1)=p1(end)/(p1(end)+p2(end)+p3(end)+p4(end));
  p(2)=p2(end)/(p1(end)+p2(end)+p3(end)+p4(end)); 
  p(3)=p3(end)/(p1(end)+p2(end)+p3(end)+p4(end));
  p(4)=p4(end)/(p1(end)+p2(end)+p3(end)+p4(end));
  rand_num=rand(1);
  criteria=0;
  if rand_num<p(1)
      emmision=emmision1;
      transmat=transmat1;
      a=find((0.2<W)&(W<0.5));
      transmat11=transmat11;
      emmision11=emmision11; 
      [emmision,transmat]=update_probabilities(all_ungap_sequences,emmision,transmat,transmat11,emmision11,a);
      emmision1=emmision;
      transmat1=transmat;
      k1=k1+1;
      [alpha1]=compute_alpha(observation,emmision1,transmat1);
      p1(k1)=sum(alpha1(:,T));
      index1=[index1,k];
      
  end
  if ((p(1)<rand_num) &&(rand_num <p(1)+p(2)))
      emmision=emmision2;
      transmat=transmat2;
      a=find((0.5<W)&(W<0.6));
      transmat11=transmat11;
      emmision11=emmision11; 
      [emmision,transmat]=update_probabilities(all_ungap_sequences,emmision,transmat,transmat11,emmision11,a);
      emmision2=emmision;
      transmat2=transmat;
      k2=k2+1;
      [alpha1]=compute_alpha(observation,emmision2,transmat2);
      p2(k2)=sum(alpha1(:,T));
      index2=[index2,k];
  end
  if  (((p(1)+p(2))<rand_num) &&(rand_num <(p(1)+p(2)+p(3))))
      emmision=emmision3;
      transmat=transmat3;
      transmat11=transmat11;
      emmision11=emmision11; 
      a=find((0.6<W)&(W<0.7));
      [emmision,transmat]=update_probabilities(all_ungap_sequences,emmision,transmat,transmat11,emmision11,a);
      emmision3=emmision;
      transmat3=transmat;
      k3=k3+1;
      [alpha1]=compute_alpha(observation,emmision3,transmat3);
      p3(k3)=sum(alpha1(:,T));
      index3=[index3,k];
  end
  if    (p(1)+p(2)+p(3))<rand_num
      emmision=emmision4;
      transmat=transmat4;
      a=find((W>0.003)&(W<0.004));
      transmat11=transmat11;
      emmision11=emmision11; 
      [emmision,transmat]=update_probabilities(all_ungap_sequences,emmision,transmat,transmat11,emmision11,a);
      emmision4=emmision;
      transmat4=transmat;
      k4=k4+1;
      [alpha1]=compute_alpha(observation,emmision4,transmat4);
      p4(k4)=sum(alpha1(:,T));
      index4=[index4,k];     
  end
end
probility_sequence=(sum(p1(find(index1>1500)))+sum(p2(find(index2>1500)))+sum(p3(find(index3>1500)))+sum(p4(find(index4>1500))))/500;
scoreml2((i-1000))=log(probility_sequence);
end