function Model_New = EM_dostep( Model, x )
%    EM_dostep.m
%      do E_step and M_step and update model parameters
%      Model = EM_dostep( Model, x )
%      input
%         Model : model structure
%         x : data to learn
%      output
%         Model : model structure updated

[ N, T ] = size( x );
m = length( Model.p );

% E-step
LL_i = zeros( m, T ); % Log Likelihood for each unit and each data
LL = zeros( T );
for i = 1:m % for each unit
  p = Model.p(i);
  x_mu  = x - repmat( p.mu, 1, T );
  x_mu2 = sum( x_mu.*x_mu, 1 );
  LL_i(i,:) = p.nu * exp( - x_mu2/p.sigma^2/2 - log( 2*pi*p.sigma^2 ) );
end
LL = sum( LL_i, 1 );

  % z_i : mean posterior for unit selection probability
  %     for each data
z_i = LL_i ./ repmat( LL, m, 1 );

% E-step update s.s.
for i = 1:m
  z = z_i( i, : );
  r.z   = 1 * sum( z ) / T;
  r.zx  = x * z' / T;
  r.zx2 = sum( x.*x, 1 ) * z' / T;
  Model.r(i) = r;
end

% M-step
for i = 1:m
  r = Model.r(i);
  mu = Model.p(i).mu;
  %Model.p(i).sigma = sqrt( (r.zx2 - 2*r.zx'*mu + r.z*mu'*mu )/N);
  Model.p(i).mu = r.zx / r.z;
end

Model_New = Model;