function [] = DST_mstep2(cur)

global TINY;

global DST;

%%LDS    
if (DST(cur).type == 3)
    
    if(cur == 1)
         S = 1;
    else
         S = DST((DST(cur).parent)).S;
    end
    
    yx_sum = zeros(DST(cur).ydim,DST(cur).xdim); 
    
    xy_sum = zeros(DST(cur).xdim,DST(cur).ydim); 
    
    yy_sum = zeros(DST(cur).ydim);
    
    P_sum = zeros(DST(cur).xdim); 
    
    A_P_sum = zeros(DST(cur).xdim,DST(cur).xdim,S);
    
    Q_P_sum = zeros(DST(cur).xdim,DST(cur).xdim,S);
    
    P2_sum = zeros(DST(cur).xdim,DST(cur).xdim,S);
    
    st_sum = zeros(S,1); 
    
    x1_sum = zeros(DST(cur).xdim,S);
    
    p1_sum = zeros(DST(cur).xdim,DST(cur).xdim,S);
    
    s1_sum = zeros(S,1);

%     Rmat = [];
%     RJ = [];
    
    for n=1:DST(1).N

        P = DST(cur).sigma{n}(:,:,1) + DST(cur).x{n}(:,1) * DST(cur).x{n}(:,1)';
        
        if(DST(cur).missing{n}(1))
            yx_sum = yx_sum + DST(cur).C * P;
            xy_sum = xy_sum + P * DST(cur).C';
            yy_sum = yy_sum + DST(cur).C * P * DST(cur).C' + DST(cur).R;
        else
            yx_sum = yx_sum + DST(cur).y{n}(1,:)' * DST(cur).x{n}(:,1)'; 
            xy_sum = xy_sum + DST(cur).x{n}(:,1) * DST(cur).y{n}(1,:);
            yy_sum = yy_sum + DST(cur).y{n}(1,:)' * DST(cur).y{n}(1,:);
        end
        
        P_sum = P_sum + P;
        
% 
%         Rmat = [Rmat ; DST(cur).y{n}(1,:) - DST(cur).x{n}(:,1)' * DST(cur).C' ; ...
%                 chol(DST(cur).sigma{n}(:,:,1) + DST(cur).x{n}(:,1) * DST(cur).x{n}(:,1)') * DST(cur).C' ; ...
%                 DST(cur).x{n}(:,1)'];
%         
%         RJ = [RJ ones(1,DST(cur).xdim+1) -1];

        for i=1:S
            
            if (cur == 1)
                s = 1;  
            else
                s = DST((DST(cur).parent)).s{n}(i,1);
            end
            
            p1_sum(:,:,i) = p1_sum(:,:,i) + s * (DST(cur).sigma{n}(:,:,1) + DST(cur).x{n}(:,1) * DST(cur).x{n}(:,1)');
            
            x1_sum(:,i) = x1_sum(:,i) + s * DST(cur).x{n}(:,1);

            s1_sum(i) = s1_sum(i) + s;

        end
        
        for t=2:DST(cur).T(n)
            
            P = DST(cur).sigma{n}(:,:,t) + DST(cur).x{n}(:,t) * DST(cur).x{n}(:,t)';
            
            if(DST(cur).missing{n}(t))
                yx_sum = yx_sum + DST(cur).C * P; 
                xy_sum = xy_sum + P * DST(cur).C';
                yy_sum = yy_sum + DST(cur).C * P * DST(cur).C' + DST(cur).R;
            else
                yx_sum = yx_sum + DST(cur).y{n}(t,:)' * DST(cur).x{n}(:,t)';     
                xy_sum = xy_sum + DST(cur).x{n}(:,t) * DST(cur).y{n}(t,:);
                yy_sum = yy_sum + DST(cur).y{n}(t,:)' * DST(cur).y{n}(t,:);
            end
       
            P_sum = P_sum + P;
            
%             Rmat = [Rmat ; DST(cur).y{n}(t,:) - DST(cur).x{n}(:,t)' * DST(cur).C' ; ...
%                 chol(DST(cur).sigma{n}(:,:,t) + DST(cur).x{n}(:,t) * DST(cur).x{n}(:,t)') * DST(cur).C' ; ...
%                 DST(cur).x{n}(:,t)'];
%         
%             RJ = [RJ ones(1,DST(cur).xdim+1) -1];
% 	
            for i=1:S
	
                if (cur == 1)
                    s = 1;  
                else
                    s = DST((DST(cur).parent)).s{n}(i,t);
                end
   
                A_P_sum(:,:,i) = A_P_sum(:,:,i) + s * (DST(cur).sigma{n}(:,:,t-1) + DST(cur).x{n}(:,t-1) * DST(cur).x{n}(:,t-1)');
            
                P2_sum(:,:,i) = P2_sum(:,:,i) + s * (DST(cur).sigma2{n}(:,:,t) + DST(cur).x{n}(:,t) * DST(cur).x{n}(:,t-1)');
                    
                Q_P_sum(:,:,i) = Q_P_sum(:,:,i) + s * (DST(cur).sigma{n}(:,:,t) + DST(cur).x{n}(:,t) * DST(cur).x{n}(:,t)');

                st_sum(i) = st_sum(i) + s;

            end
            
        end
    
    end
     
      DST(cur).C = yx_sum * safeinv(P_sum);
    

%%%Good    
 
      DST(cur).R = ( yy_sum - (DST(cur).C * xy_sum) - (DST(cur).C * xy_sum)' + (DST(cur).C * P_sum * DST(cur).C')) / (sum(DST(cur).T));

      for i=1:S
 
         DST(cur).A(:,:,i) = P2_sum(:,:,i) * safeinv(A_P_sum(:,:,i) ); 
              
        %%Good
        
         DST(cur).Q(:,:,i) = (Q_P_sum(:,:,i) - DST(cur).A(:,:,i) * P2_sum(:,:,i)' - P2_sum(:,:,i) * DST(cur).A(:,:,i)'...
             + DST(cur).A(:,:,i) * A_P_sum(:,:,i) * DST(cur).A(:,:,i)') / st_sum(i)  + 1e-8 * eye(DST(cur).xdim);
        
        
         DST(cur).Qinv(:,:,i) = safeinv(DST(cur).Q(:,:,i)); 

         DST(cur).mu(:,i) = x1_sum(:,i)/s1_sum(i);
          
        %%%Good
         DST(cur).Q1(:,:,i) = (p1_sum(:,:,i) - x1_sum(:,i) * DST(cur).mu(:,i)' - DST(cur).mu(:,i) * x1_sum(:,i)' )/ s1_sum(i) + ...
             DST(cur).mu(:,i) * DST(cur).mu(:,i)' + 1e-8 * eye(DST(cur).xdim);
        

         DST(cur).Q1inv(:,:,i) = safeinv(DST(cur).Q1(:,:,i)); 

     end
     
%%HMM or internal chain
else 

    if (cur==1)
        pS = 1;
    else
        pS = DST(DST(cur).parent).S;
   
    end
    
        DST(cur).Phi = zeros(DST(cur).S,DST(cur).S,pS);
        
        DST(cur).phi = zeros(DST(cur).S,pS);
        
        for k=1:pS
      
            for n=1:DST(1).N
                
                if (cur==1)
                    ps = ones(1,DST(cur).T(n));
                else
                    ps = DST(DST(cur).parent).s{n};
                end
		
                for t=2:DST(cur).T(n)
                    

                    DST(cur).Phi(:,:,k) = DST(cur).Phi(:,:,k) + DST(cur).xi{n}(:,:,t-1) * ps(k,t);
		
                end
                
                DST(cur).phi(:,k) = DST(cur).phi(:,k) + DST(cur).s{n}(:,1) * ps(k,1);

            end
            
            DST(cur).phi(:,k) = normalise2(DST(cur).phi(:,k),1e-25);
            
            DST(cur).Phi(:,:,k) = normalise2(DST(cur).Phi(:,:,k),1e-25);

%                		
%             DST(cur).phi(:,k) = normalise( DST(cur).phi(:,k));
%             
%             DST(cur).Phi(:,:,k) = mk_stochastic(DST(cur).Phi(:,:,k)')';
            
        end

    %%HMM learn emission parameters
    %fix sometime
    if(DST(cur).type==2)
        
        DST(cur).eta = zeros(DST(cur).S,DST(cur).S);
        
        for t=1:T
            o = DST(cur).y(t);
         
            DST(cur).eta(o,:) = DST(cur).eta(o,:) + DST(cur).s(:,t)';
        end
    
       DST(cur).eta = mk_stochastic(DST(cur).eta')';
            
       
    %else internal node, recurse
    else
        
        for m=1:length(DST(cur).child)
            
            DST_mstep2(DST(cur).child(m));
            
        end
        
    end
       
end