function [param] = em(X,param);

[n,m] = size(X);
k = length(param);

logp = zeros(n,k);
for i=1:k,
 logp(:,i) = evalgauss(X,param(i));
end;

pt = softmax(logp); % posterior component assignments 

for i=1:k,
  ni = sum(pt(:,i));
  param(i).p = ni/n;
  param(i).mean = sum( (pt(:,i)*ones(1,m)).*X )/ni; % weighted mean
  Xt = (X-ones(n,1)*param(i).mean);
  param(i).cov = ((pt(:,i)*ones(1,m)).*Xt)'*Xt/ni; % weighted covariance
end;



