gmm_test.m

来自「patten regnization source从1-14章能运行」· M 代码 · 共 47 行

% 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;