Gaussian mixture model (MATLAB code + comment)

Here I am studying the EMGMM code in statistical patte7rn recognition Toolbox, where the main points of knowledge in the code are explained in the previous GMM document, including the two steps in the EM algorithm. I first look at the principle, and then to see the code, in the process of commenting on the code I re-combed the whole theory system, or is very grateful to this way, do one thing.

Main function emgmm
Child functions
Probability density estimation of pdfgauss multivariate Gaussian distribution
Data classification method used during Knnrule Knnclass model parameter initialization
MLCGMM estimation of GMM parameters known to data classes
M-Step parameter estimation in MELGMM EM algorithm

function model=emgmm (x,options,init_model)% emgmm expectation-maximization algorithm for Gaussian mixture model.  % Synopsis:% model = EMGMM (X)% model = EMGMM (x,options)% model = EMGMM (x,options,init_model)% Description:% This function implements the Expectation-maximization algorithm% (EM) [schles68][dlr77] which computes the maximum-like 
Lihood% estimate of the paramaters of the Gaussian mixture model (GMM). % The EM algorithm is a iterative procedure which monotonically% increases log-likelihood of the current estimate unt 
Il it reaches% a local optimum.
% of the number of components of the GMM are given in options.ncomp% (default 2). %%%%%em algorithm iteration stop condition% The following three stopping is condition used:% 1. Improvement of the Log-likelihood is less than given% threshold% LOGL (t+1)-LOGL (t) < OPTIONS.E PS_LOGL%%%LOGL (t) is monotonically increasing, and as the lower bound increases, it gradually approaches its maximum value.

Stop iterating when it no longer changes or changes in the range of hours. % 2. Change of the squared differences of a estimated Posteriory% probabilities is less than given threshold% | | Alpha (T+1)-alpha (t) | | ^2 < Options.eps_alpha%%%%e the desired of the implied variable by the given parameter, m step to re-estimate the parameters according to expectation, return to e-step, and draw the expectation of the new implied variable.

By finding the L2 norm of both vectors, the iteration is stopped when its%%%% is constant or changes very small. % 3.
Number of iterations exceeds given threshold.    % T >= Options.tmax%%%% set maximum iteration count% The type of estimated covariance matrices is optional:% Options.cov_type = ' full ' full covariance matrix (default)% Options.cov_type = ' diag ' diagonal covarinace m Atrix% Cov_options.type = ' spherical ' spherical covariance matrix%%%% in model parameter initialization, KNN is used to classify the sample data and the maximum likelihood estimation method is used to estimate the parameters, the following three types are KNN Center Point initialization method% The initial model (estimate) is selected:% 1. Randomly (options.init = ' random ')% 2. Using C-means (options.init = ' Cmeans ')% 3.
Using the user specified Init_model.
% Input:% x [Dim x Num_data] Data sample.
% options [struct] Control paramaters:%. Ncomp [1x1] Number of components of GMM (default 2). %  . Tmax [1x1] maximal number of iterations (default INF).
%. EPS_LOGL [1x1] Minimal improvement in Log-likelihood (default 0).
%. Eps_alpha [1x1] Minimal change of Alphas (default 0).
%. Cov_type [1x1] type of estimated covarince matrices (see above).
%. init [string] ' random ' use random initial model (default);
% ' Cmeans ' use K-means to find initial model.
%. verb [1x1] If 1 then info is displayed (default 0). % Init_model [struct] Initial model:%.
Mean [Dim x Ncomp] Mean vectors. %   .
Cov [Dim x Dim x Ncomp] Covariance matrices. %   .
Priors [1 x ncomp] Weights of mixture components. %   .
Alpha [Ncomp x Num_data] (optional) distribution of hidden state.
%. t [1x1] (optional) Counter of iterations. % Output:% model [struct] estimated Gaussian mixture model:%.
Mean [Dim x Ncomp] Mean vectors. %   .
Cov [Dim x Dim x Ncomp] Covariance matrices. %   .
Prior [1 x ncomp] Weights of mixture components.
%. t [1x1] number iterations. %. options [StruCT] Copy of used options.
%. exitflag [int] 0 ... maximal number of iterations was exceeded. % 1 or 2 ... EM has converged;
Indicates which stopping% was used (see above).
% Example:% note:if EM algorithm does not converge run it again from different% initial model. % EM is used to estimate parameters of mixture of 2 Guassians:% True_model = struct (' Mean ', [-2 2], ' Cov ', [1 0.5], ' Prior
', [0.4 0.6]);
% Sample = Gmmsamp (True_model, 100); % Estimated_model = EMGMM (sample.
X,struct (' Ncomp ', 2, ' verb ', 1)); %% Figure; Ppatterns (Sample.
X);
% h1=pgmm (true_model,struct (' Color ', ' r '));
% h2=pgmm (estimated_model,struct (' Color ', ' B ')); 
% Legend ([H1 (1) H2 (1)], ' Ground truth ', ' ML estimation '); % Figure; Hold on; Xlabel (' iterations ');
Ylabel (' Log-likelihood ');
% plot (ESTIMATED_MODEL.LOGL);
% See also% mlcgmm, Mmgauss, PDFGMM, Gmmsamp. % about:statistical patte7rn recognition Toolbox% (C) 1999-2003, written by Vojtech FranC and Vaclav Hlavac% <a href= "http://www.cvut.cz" >czech Technical University prague</a>% <a href= "/http www.feld.cvut.cz ">faculty of electrical engineering</a>% <a href=" http://cmp.felk.cvut.cz ">center for Machine perception</a>% Modifications:% 26-may-2004, VF, initialization by K-means added% 1-may-2004, VF% 19-se  p-2003, VF% 16-mar-2003, VF% processing input arguments%-----------------------------------------if Nargin < 2, Options=[]; else Options=c2s (options); End If ~isfield (options, ' Ncomp '), Options.ncomp = 2; End If ~isfield (options, ' Tmax '), Options.tmax = 1000; End If ~isfield (options, ' Eps_alpha '), Options.eps_alpha = 0; End If ~isfield (options, ' Eps_logl '), OPTIONS.EPS_LOGL = 0; End If ~isfield (options, ' Cov_type '), Options.cov_type = ' full '; End If ~isfield (options, ' init '), Options.init = ' Cmeans '; End If ~isfield (options, ' verb '), options.verb = 0;      End [Dim,num_data] = size (X); %%%%% data input format to be aware of the% setupInitial model%---------------------------------if Nargin = = 3,% take model from input%------------------------ 
  -----model = Init_model; If ~isfield (model, ' t '), model.t = 0; End If ~isfield (model, ' Alpha '), model.
  Alpha=-inf*ones (Options.num_gauss,num_data); End If ~isfield (model, ' LOGL '), Model.logl=-inf;
    End else% compute initial model%------------------------------------switch Options.init,% random model Case ' random '% takes randomly first Num_gauss trn.  
     Vectors as mean vectors inx = randperm (Num_data);
     Inx=inx (1:options.ncomp);

    Centers_x = X (:, Inx);
     % K-means Clustering case ' cmeans ' tmp = Cmeans (X, Options.ncomp); centers_x = tmp.

    X
  Otherwise error (' Unknown initialization method. ');
  End KNN = Knnrule ({' X ', centers_x, ' Y ', [1:options.ncomp]},1);                                         y = Knnclass (X,KNN); %%%%%% author default 1 nearest neighbor classifies data into n categories (n component)% uses ML estimation of complete DATA model = MLCGMM ({' X ', x, ' Y ', y}, Options.cov_type); The%%%%%% data category is known to use the maximum likelihood estimation model.
  Alpha = zeros (options.ncomp,num_data); For i = 1:options.ncomp, model.                               Alpha (I,find (y==i)) = 1;
  %%%% marks the K line of the corresponding category k in each column as 1 end model.logl=-inf;
  model.t = 1;
  model.options = options;

Model.fun = ' pdfgmm ';
End% Main loop of EM algorithm%-------------------------------------model.exitflag = 0;

  While Model.exitflag = = 0 & model.t < Options.tmax,% counter of iterations model.t = model.t + 1;
  %----------------------------------------------------% e-step The distribution of hidden states is computed based
  % on the current estimate.
  %----------------------------------------------------A=pdfgauss (X, model); Newalpha = (model.          Prior (:) *ones (1,num_data)). *pdfgauss (X, model);                                       %%%% The estimated value NEWLOGL = SUM (log (sum (newalpha,1)) according to the E-step expectation formula; %%%%% to find the new likelihood probability value Newalpha = newalpha./(Ones (options.ncomp,1) *sum (newalpha,1));
  %%%%% the sum of the probabilities of satisfying data points in each conponent is 1 b=a.*newalpha;
  %------------------------------------------------------% stopping conditions. %------------------------------------------------------% 1) change in distribution of hidden state Alpha Model.delta _alpha = SUM (SUM (model.                 Alpha-newalpha). ^2));                              %%%% update hidden state distribution% 2) change in Log-likelihood MODEL.DELTA_LOGL = NEWLOGL-MODEL.LOGL (end);

  %%%% the new calculated likelihood probability value minus the last likelihood probability estimate, guaranteeing MODEL.LOGL = [Model.logl NEWLOGL] within the iteration termination condition; If Options.verb, fprintf ('%d:logl=%f, delta_logl=%f, delta_alpha=%f\n ',... model.t, MODEL.LOGL (end), model.de
  LTA_LOGL, Model.delta_alpha);
  End If Options.eps_logl >= model.delta_logl, model.exitflag = 1;
  ElseIf options.eps_alpha >= model.delta_alpha, model.exitflag = 2; else model.

    Alpha = Newalpha; %----------------------------------------------------% m-step the new parameters maximizing expectation of% Log-likelihood are computed. %----------------------------------------------------Tmp_model = melgmm (X,model.                       Alpha,options.cov_type); %%%%% based on the new implicit state distribution estimation parameters, the same maximum likelihood estimation method is used to calculate the estimate model according to the formula in M step. Mean = Tmp_model.
    Mean; Model. Cov = Tmp_model.
    Cov; Model. Prior = Tmp_model.

  Prior;

End end% while main loop return;
 NOTE: Child functions are no longer posted, if you need a complete set of code, you can download the Toolbox http://cmp.felk.cvut.cz/cmp/software/stprtool/