function compParams = Mstep(X, p)
  [n,d] = size(X) ;
  k = size(p,2) ;
  compParams = [] ;
  for i=1:k
    reped_p = repmat(p(:,i),1,d) ;	% This can come in usefull
    ni = sum(p(:,i));			% Total 'weight' of ith component
    prior = ni/n;
    mu = sum( reped_p.*X )/ni; % weighted mean
    XminusMu = op(X,'-',mu);
    Sigma = (reped_p.*XminusMu)'*XminusMu/ni; % weighted covariance
    compParams = [compParams, struct('prior',prior,'mean',mu,'cov',Sigma)];
  end;

    