% EM 舅絆府硫欄肺 彌利狼 去欽 啊快矯救 葛膽傅 
% 券屈 單撈磐 肺爹竅扁 
load train1; 
% 檬扁 頗扼固磐 汲瀝 
[n D] = size(x);     % n : 包螟 單撈磐 薦, D : 瞞盔
k = 6;               % 去欽薦 
p = ones(1,k)/k;     % 去欽 厚啦 檬扁藹
mu = randn(D,k);     % 乞閉 
s2 = zeros(D,D,k);   % 傍盒魂 青紡
niter=100;           % 館汗 冉薦 
% 措阿 己盒惑俊輯 瘤薦利欄肺 刀賦茄 盒魂藹 檬扁拳 
for i=1:k
  s2(:,:,i) = -100*diag(log(rand(D,1))); % variances
end

set(gcf,'Renderer','zbuffer');

clear Z;
try
  % niter 館汗竅咯 EM 切嚼  
  for t=1:niter,
    fprintf('t=%d\r',t);
    % E-竄拌:
    for i=1:k
      Z(:,i) = p(i)*det(s2(:,:,i))^(-0.5)*exp(-0.5*sum((x'-repmat(mu(:,i),1,n))'*inv(s2(:,:,i)).*(x'-repmat(mu(:,i),1,n))',2));
    end
    Z = Z./repmat(sum(Z,2),1,k);
    
    % M-竄拌:
    for i=1:k
      mu(:,i) = (x'*Z(:,i))./sum(Z(:,i));
      s2(:,:,i) = (x'-repmat(mu(:,i),1,n))*(repmat(Z(:,i),1,D).*(x'-repmat(mu(:,i),1,n))')./sum(Z(:,i));
      p(i) = mean(Z(:,i));
    end
    
    clf
    hold on
    plot3(x(:,1),x(:,2),x(:,3),'.');
    for i=1:k
      plot_gaussian(s2(:,:,i),mu(:,i),i,20);
    end
    drawnow;
  end
catch
  disp('風橇惑俊輯 薦摹利 坷幅 慣積 - 酒付檔 漂撈 青紡撈 慣積茄 巴 鞍瀾');
end;