traingmm.m

patten regnization source从1-14章能运行

% [mus, sigmas, weights] = trainGMM(X, mus, sigmas, weights, threshold, maxiterations, modeltype);
% 单捞磐 青纺 X惑俊辑 [mus, sigmas, weights]牢 颇扼固磐俊 狼窍咯 沥秦柳 GMM阑  
% EM 舅绊府硫阑 捞侩茄 切嚼阑 肺弊快档啊 THRESHOLD 蔼焊促 利芭唱 
% MAXITERATIONS 馆汗 寝荐啊 瞪 锭 鳖瘤 荐青.
% 去钦 己盒狼 傍盒魂俊 荤侩等 葛胆狼 屈怕绰 MODELTYPE俊辑 
% 胶纳老等 窜困 青纺 趣篮 措阿 青纺肺 瘤沥 啊瓷. 
% (肯傈 傍盒魂栏肺狼 犬厘且妨搁 sigmas啊 氦磐啊 酒聪扼 青纺肺 沥秦廉具父 茄促.) 
%
% 涝仿 牢磊 :
% X绰 青捞 d瞒盔牢 漂隆 氦磐 单捞磐 青纺.
% mus : d*M 青纺肺 弊 青篮 雇胶媚 己盒俊 措茄 乞闭 氦磐.
% sigmas : d*M 青纺肺 弊 青篮 雇胶媚 去钦 己盒俊 措茄 盒魂 氦磐.
% weights : M-瞒盔 氦磐肺 去钦 啊吝摹.
% modeltype : 'scalarcov' 趣篮 'diagcov' 巩磊凯.
% 免仿 牢磊 : 

function [mus, sigmas, weights] = trainGMM(X, mus, sigmas, weights, threshold, maxiterations, modeltype);
  
% Check input params and get dimensions:
  [d,numpoints] = size(X);
  M = length(weights);
  if (abs(sum(weights) - 1) < 0.00001) == 0
     error('trainGMM: Weights should sum to 1!');
  end
  
  if (size(mus) == [d,M] & size(sigmas) == [d,M]) == 0
     error('trainGMM: Weights should sum to 1!');
  end
 
  if strcmp(modeltype,'scalarcov')
    xvars = zeros(1,numpoints);
  elseif strcmp(modeltype,'diagcov')
    xvarmat = zeros(d,numpoints);
  else
    error('Unknown model type!');
  end

% Plot data:
  % hold off
  % plot(X(1,:),X(2,:),'r.')
  % hold on

% Plot initial component centroids:
  %  plot(mus(1,:),mus(2,:),'bo')

% Calculate initial likelihoods:
  [loglikevec, pmat] = loglikeGMM(X, mus, sigmas, weights);
  loglike = sum(loglikevec);
  disp('Log likelihood at initialisation: '), disp(loglike);

% Iterate until convergence:
  for iteration = 1:maxiterations
    disp('Iteration:'),disp(iteration);

% E step: calculate posteriors under current parameter values:
    for n = 1:numpoints
      pmat(:,n) = pmat(:,n)/sum(pmat(:,n));
    end
    disp(pmat(:,1:5))

% M step: update parameter values to maximise posteriors:
    for m = 1:M;
% Mean update:
      psum = sum(pmat(m,:));
      mus(:,m) = (X * pmat(m,:)')/psum;
% Variance update:
      if strcmp(modeltype,'scalarcov')
	for n = 1:numpoints
          meansub = X(:,n)-mus(:,m);
          xvars(n) = meansub'*meansub;
	end
	sigmas(:,m) = ones(d,1)*(xvars * pmat(m,:)')/(d*psum);
      elseif strcmp(modeltype,'diagcov')
	for n = 1:numpoints
          meansub = X(:,n)-mus(:,m);
          xvarmat(:,n) = meansub.*meansub;
	end
	sigmas(:,m) = xvarmat * pmat(m,:)' / psum;       
      end
% Weight update:
      weights(m) = psum/numpoints;
    end

% Plot new component centroids:
  drawnow
  plot(mus(1,:),mus(2,:),'bo') % 颇扼固磐 弥利拳 苞沥阑 焊咯林扁 困秦辑
  disp(weights)
  disp(mus)
  disp(sigmas)

% Calculate new log likelihood:
    [newloglikevec, pmat] = loglikeGMM(X, mus, sigmas, weights);
    newloglike = sum(newloglikevec);
    disp('Log likelihood: '), disp(newloglike);

% Test for convergence:
    if (newloglike < loglike)
      error('Log likelihood decreased!')
    end
    if ((newloglike - loglike)/numpoints < threshold)
      disp('Converged.')
      break;
    else 
      loglike = newloglike;
    end;

  end
  disp(' '),disp(' '),disp(' ')