Download All

Code covered by the BSD License  

Highlights from
Statistical Learning Toolbox

from Statistical Learning Toolbox by Dahua Lin
Functions for statistical learning, pattern recognition and computer vision, covering many topics.

Description of slfmm
Home > sltoolbox > stat > slfmm.m

slfmm

PURPOSE ^

SLFMM Learns a Finite Mixture Model (FMM)

SYNOPSIS ^

function [S, cw, pp, info] = slfmm(X, n, estfunctor, evalfunctor, varargin)

DESCRIPTION ^

SLFMM Learns a Finite Mixture Model (FMM)

 $ Syntax $
   - [S, cw] = slfmm(X, n, estfunctor, evalfunctor, ...)
   - [S, cw, pp] = slfmm(X, n, estfunctor, evalfunctor, ...)
   - [S, cw, pp, info] = slfmm(X, n, estfunctor, evalfunctor, ...)
 
 $ Arguments $
   - X:            the samples
   - n:            the number of samples in X
   - estfunctor:   the functor to estimate the model from weighted 
                   samples, it should take the following form.
                   The functor is a cell array, with the first element
                   being a function, while the others are the extra
                   parameters: {estfunc, ...}
                   For different modes, the function should take different
                   forms:
                   - 'simple' mode:
                       s = estfunc(s0, X, n, w, ...)
                   - 'innermul' mode:
                       S = estfunc(S0, X, n, W, selinds, ...)
                   Here
                   - s:    a single model
                   - s0:   the previous single model (initially, it is [])
                   - S:    a multi-model
                   - S0:   the previous multi model (initially, it is [])
                   - X:    the sample set
                   - n:    the number of samples
                   - w:    the 1 x n weight vector for a specific
                           component
                   - W:    the k x n weight matrix for all components
                   - selinds:  the selected indices of models to be
                               updated. If it is empty, all models need
                               to be updated.
   - evalfunctor:  the functor to evaluate component-conditional 
                   likelihood
                   condp = evalfunc(S, X, n, ...)
                   for 'simple' mode, condp is a length-n vector
                   for 'innermul', condp is k x n matrix.
                   If condpmode is 'log', then the condp is logarithm of
                   the probabilities.  
   - S:        the struct of the learned finite mixtured model
   - cw:       the component mixture weights
   - pp:       the posteriori of each sample w.r.t mixture components

 $ Description $
   - [S, cw, pp] = slfmm(X, n, estfunc, evalfunc, ...) learns the finite 
     mixture model from the samples in X, according to the properties.
     \*
     \t  Table 1. Properties of GMM Learning                 \\
     \h    name      &      description                      \\
          'method'   & The method using for GMM learning. Currently,
                       there are only one method available: 'EM',
                       default = 'EM'.
          'update'   & The way of updating (default = 'pass'):
                        1. 'pass'     Pass-wise update;
                        2. 'comp'     Component-wise update     \\
          'cyclecn'  & The ratio of components to be updated in each 
                       cycle for 'comp' update scheme. default = 1.
          'iter'     & The iteration control properties for sliterproc
                       default = {'maxiter', 100}
          'tol'      & The maximum tolerance of posteriori error when 
                       the iteration terminates, (default = 1e-6) \\
          'verbose'  & whether to display information while iteration, 
                       default = true                              \\
          'initc'    & The labels of initial clustering, if not specified, 
                       we will use random clustering. If initc is specified,
                       K will be forced to the the number of unique labels.
                       default = []. \\
          'initpp'   & The initial posteriori of samples. (k0 x n)  
                       default = []. \\
          'annthres' & The threshold of annealing in FJ algorithm, 
                       default = 0 (unit = average mixture weight) \\
          'weights'  & The weights of the samples, default = [], 
                       indicating non-weighted. Weights should be given
                       by a 1 x n row vector. \\
          'estmode'  & The mode of estimation:
                       'simple':    S is a struct array of models, each
                                    element represents a model, there is
                                    no internal representation of 
                                    component weights. (default)
                       'innermul':  S is a struct or an object maintaining
                                    a set of models
          'condpmode' & The mode of evaluated likelihood
                        'normal':  they are the likelihood values
                        'log':     they are the logarithm of the
                                   likelihood
          'manifunc'  & The function to manipulate the models
                        (only for innermul mode)
                        S = manifunc(S0, 'select', [k1, k2, ...])
                           to retain the k1-th, k2-th, ... models
     \*

   - [S, cw, pp, info] = slfmm(X, n, estfunc, evalfunc, ...)
     also outputs the information of learning process. The info is a 
     struct with following fields
       - numiters:     the number of iterations
       - stopdiscrep:  the discrepance value on stop
       - isconverged:  whether the process has been converged

 $ Remarks $
   - Note that there is no any restriction on the form of X. The
     only condition is that the estfunc and evalfunc accept it.

   - When annthres is larger than zero, the manifunc supporing
     the operation 'discard' is required.

   - For initialization, you should specify either of the initc or 
     initpp. If both initc or initpp are given, only the initpp
     will take effect.
   
 $ History $
   - Created by Dahua Lin on Aug 17, 2006
       - as a foundation of all finite mixture models, such as GMM
       - it is based on the original slgmm function
   - Modified by Dahua Lin on Aug 29, 2006
       - based on the sliterproc to control iteration process.
   - Modified by Dahua Lin on Sep 10, 2006
       - replace slmul by slmulvec to increase efficiency

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:
  • slgmm SLGMM Learns Gaussian Mixture model from samples

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [S, cw, pp, info] = slfmm(X, n, estfunctor, evalfunctor, varargin)
0002 %SLFMM Learns a Finite Mixture Model (FMM)
0003 %
0004 % $ Syntax $
0005 %   - [S, cw] = slfmm(X, n, estfunctor, evalfunctor, ...)
0006 %   - [S, cw, pp] = slfmm(X, n, estfunctor, evalfunctor, ...)
0007 %   - [S, cw, pp, info] = slfmm(X, n, estfunctor, evalfunctor, ...)
0008 %
0009 % $ Arguments $
0010 %   - X:            the samples
0011 %   - n:            the number of samples in X
0012 %   - estfunctor:   the functor to estimate the model from weighted
0013 %                   samples, it should take the following form.
0014 %                   The functor is a cell array, with the first element
0015 %                   being a function, while the others are the extra
0016 %                   parameters: {estfunc, ...}
0017 %                   For different modes, the function should take different
0018 %                   forms:
0019 %                   - 'simple' mode:
0020 %                       s = estfunc(s0, X, n, w, ...)
0021 %                   - 'innermul' mode:
0022 %                       S = estfunc(S0, X, n, W, selinds, ...)
0023 %                   Here
0024 %                   - s:    a single model
0025 %                   - s0:   the previous single model (initially, it is [])
0026 %                   - S:    a multi-model
0027 %                   - S0:   the previous multi model (initially, it is [])
0028 %                   - X:    the sample set
0029 %                   - n:    the number of samples
0030 %                   - w:    the 1 x n weight vector for a specific
0031 %                           component
0032 %                   - W:    the k x n weight matrix for all components
0033 %                   - selinds:  the selected indices of models to be
0034 %                               updated. If it is empty, all models need
0035 %                               to be updated.
0036 %   - evalfunctor:  the functor to evaluate component-conditional
0037 %                   likelihood
0038 %                   condp = evalfunc(S, X, n, ...)
0039 %                   for 'simple' mode, condp is a length-n vector
0040 %                   for 'innermul', condp is k x n matrix.
0041 %                   If condpmode is 'log', then the condp is logarithm of
0042 %                   the probabilities.
0043 %   - S:        the struct of the learned finite mixtured model
0044 %   - cw:       the component mixture weights
0045 %   - pp:       the posteriori of each sample w.r.t mixture components
0046 %
0047 % $ Description $
0048 %   - [S, cw, pp] = slfmm(X, n, estfunc, evalfunc, ...) learns the finite
0049 %     mixture model from the samples in X, according to the properties.
0050 %     \*
0051 %     \t  Table 1. Properties of GMM Learning                 \\
0052 %     \h    name      &      description                      \\
0053 %          'method'   & The method using for GMM learning. Currently,
0054 %                       there are only one method available: 'EM',
0055 %                       default = 'EM'.
0056 %          'update'   & The way of updating (default = 'pass'):
0057 %                        1. 'pass'     Pass-wise update;
0058 %                        2. 'comp'     Component-wise update     \\
0059 %          'cyclecn'  & The ratio of components to be updated in each
0060 %                       cycle for 'comp' update scheme. default = 1.
0061 %          'iter'     & The iteration control properties for sliterproc
0062 %                       default = {'maxiter', 100}
0063 %          'tol'      & The maximum tolerance of posteriori error when
0064 %                       the iteration terminates, (default = 1e-6) \\
0065 %          'verbose'  & whether to display information while iteration,
0066 %                       default = true                              \\
0067 %          'initc'    & The labels of initial clustering, if not specified,
0068 %                       we will use random clustering. If initc is specified,
0069 %                       K will be forced to the the number of unique labels.
0070 %                       default = []. \\
0071 %          'initpp'   & The initial posteriori of samples. (k0 x n)
0072 %                       default = []. \\
0073 %          'annthres' & The threshold of annealing in FJ algorithm,
0074 %                       default = 0 (unit = average mixture weight) \\
0075 %          'weights'  & The weights of the samples, default = [],
0076 %                       indicating non-weighted. Weights should be given
0077 %                       by a 1 x n row vector. \\
0078 %          'estmode'  & The mode of estimation:
0079 %                       'simple':    S is a struct array of models, each
0080 %                                    element represents a model, there is
0081 %                                    no internal representation of
0082 %                                    component weights. (default)
0083 %                       'innermul':  S is a struct or an object maintaining
0084 %                                    a set of models
0085 %          'condpmode' & The mode of evaluated likelihood
0086 %                        'normal':  they are the likelihood values
0087 %                        'log':     they are the logarithm of the
0088 %                                   likelihood
0089 %          'manifunc'  & The function to manipulate the models
0090 %                        (only for innermul mode)
0091 %                        S = manifunc(S0, 'select', [k1, k2, ...])
0092 %                           to retain the k1-th, k2-th, ... models
0093 %     \*
0094 %
0095 %   - [S, cw, pp, info] = slfmm(X, n, estfunc, evalfunc, ...)
0096 %     also outputs the information of learning process. The info is a
0097 %     struct with following fields
0098 %       - numiters:     the number of iterations
0099 %       - stopdiscrep:  the discrepance value on stop
0100 %       - isconverged:  whether the process has been converged
0101 %
0102 % $ Remarks $
0103 %   - Note that there is no any restriction on the form of X. The
0104 %     only condition is that the estfunc and evalfunc accept it.
0105 %
0106 %   - When annthres is larger than zero, the manifunc supporing
0107 %     the operation 'discard' is required.
0108 %
0109 %   - For initialization, you should specify either of the initc or
0110 %     initpp. If both initc or initpp are given, only the initpp
0111 %     will take effect.
0112 %
0113 % $ History $
0114 %   - Created by Dahua Lin on Aug 17, 2006
0115 %       - as a foundation of all finite mixture models, such as GMM
0116 %       - it is based on the original slgmm function
0117 %   - Modified by Dahua Lin on Aug 29, 2006
0118 %       - based on the sliterproc to control iteration process.
0119 %   - Modified by Dahua Lin on Sep 10, 2006
0120 %       - replace slmul by slmulvec to increase efficiency
0121 %
0122 
0123 %% parse and verify input arguments
0124 
0125 props.method = 'EM';
0126 props.update = 'pass';
0127 props.cyclecn = 1;
0128 props.iter = {'maxiter', 100};
0129 props.tol = 1e-6;
0130 props.verbose = true;
0131 props.initc = [];
0132 props.initpp = [];
0133 props.annthres = 0;
0134 props.weights = [];
0135 props.estmode = 'simple';
0136 props.condpmode = 'normal';
0137 props.manifunc = [];
0138 
0139 props = slparseprops(props, varargin{:});
0140 
0141 % add two properties
0142 props.estfunctor = estfunctor;
0143 props.evalfunctor = evalfunctor;
0144 
0145 
0146 checkvalid('learning method', props.method, {'EM'});
0147 checkvalid('updating scheme', props.update, {'pass', 'comp'});
0148 checkvalid('estimation mode', props.estmode, {'simple', 'innermul'});
0149 checkvalid('conditional probability form', props.condpmode, {'normal', 'log'});
0150 
0151 if isempty(props.initc) && isempty(props.initpp)
0152     error('sltoolbox:invalidarg', ...
0153         'Please specify either initc or initpp');
0154 end
0155 
0156 if ~isempty(props.weights)
0157     if ~isequal(size(props.weights), [1, n])
0158         error('sltoolbox:sizmismatch', ...
0159             'The weights should be a 1 x n row vector');
0160     end
0161 end
0162 
0163 if props.annthres > 0 
0164     if isempty(props.manifunc)
0165         error('sltoolbox:rterror', ...
0166             'The manipulation function is required when component annealing is on');
0167     end
0168 end
0169 
0170 
0171 %% initialization
0172 
0173 slsharedisp_attach('slfmm', 'show', props.verbose);
0174 
0175 % determine initial K value and initial clustering
0176 if ~isempty(props.initpp)
0177     
0178     pp = props.initpp;
0179     W = weightmap(pp, props.weights);
0180     
0181 elseif ~isempty(props.initc)
0182     
0183     [pp, W] = labels2weights(props.initc, props.weights);
0184         
0185 end
0186 
0187 initK = size(W, 1);
0188 
0189 slsharedisp('FMM Learning Parameters:');
0190 slsharedisp_incindent;
0191 slsharedisp('method  = %s', props.method);
0192 slsharedisp('update  = %s', props.update);
0193 slsharedisp('estmode = %s', props.estmode);
0194 slsharedisp('init K  = %d', initK);
0195 slsharedisp('anneal thres = %g', props.annthres);
0196 slsharedisp_decindent;
0197 
0198 %% updating
0199 
0200 switch props.method       
0201     case 'EM'
0202         slsharedisp('Run Finite-Mixture Model Learning by EM');
0203         
0204         em_estfunctor = {@fmm_em_est, props};
0205         em_evalfunctor = {@fmm_em_eval, props};
0206         em_cmpfunctor = {@fmm_em_cmp, props};
0207         
0208         models = {[], []};
0209         data = {X, n};
0210         Q = {pp, W};
0211         
0212         props.iter = {props.iter{:}, 'titlebreak', false};
0213         [models, Q, info] = slreevallearn(models, Q, data, ...
0214             em_estfunctor, em_evalfunctor, em_cmpfunctor, ...
0215             'iter', props.iter, 'isrecorded', false);
0216         
0217         S = models{1};
0218         cw = models{2};
0219         pp = Q{1};
0220     
0221 end
0222 
0223 slsharedisp_detach;
0224 
0225 
0226 %% The Iteration functions for EM
0227 
0228 % data = {X, n}
0229 % Q = {pp, W}
0230 % models = {S, cw}
0231 
0232 function models = fmm_em_est(models, data, Q, props)
0233 
0234 S = models{1};
0235 X = data{1};  
0236 % n = data{2};
0237 W = Q{2};
0238 
0239 if isempty(S)  % initialize model
0240     [S, cw] = init_models(X, W, props);
0241     
0242 else           % update model
0243     switch props.update
0244         case 'pass'
0245             [S, W, cw] = update_models(S, X, W, [], props);
0246             slignorevars(W);
0247 
0248         case 'comp'
0249             cn = max(ceil(props.cyclecn * size(W, 1)), 1);
0250             for t = 1 : cn
0251                 si = ceil(rand * size(W,1));
0252                 si = min(max(si, 1), size(W, 1));
0253                 [S, W, cw] = update_models(S, X, W, si, props);
0254             end
0255     end
0256 end
0257 
0258 models = {S, cw};
0259 
0260 
0261 function Q = fmm_em_eval(models, data, Q, props)
0262 
0263 S = models{1}; cw = models{2};
0264 X = data{1}; n = data{2};
0265 slignorevars(Q);
0266 
0267 [pp, W] = update_weightmap(S, cw, X, n, props);
0268 
0269 Q = {pp, W};
0270 
0271 
0272 function isconverged = fmm_em_cmp(models_prev, models, Q_prev, Q, props)
0273 
0274 slignorevars(models_prev, models);
0275 
0276 W_prev = Q_prev{2};
0277 W = Q{2};
0278 isconverged = false;
0279 
0280 slsharedisp_attach('fmm_em_cmp');
0281 
0282 if isequal(size(W), size(W_prev))
0283     wdiff = sldiscrep(W_prev, W, 'maxdiffnrm', true);
0284     slsharedisp('K = %d: wdiff = %g', size(W, 1), wdiff);
0285     if wdiff < props.tol
0286         isconverged = true;
0287     end
0288 else
0289     slsharedisp('K = %d -> %d', size(W_prev, 1), size(W, 1));
0290 end
0291 
0292 slsharedisp_detach;
0293 
0294 
0295 
0296 
0297 %% Core routines
0298 
0299 % The function to initialize models
0300 function [S, cw] = init_models(X, W, props)
0301 
0302 [k, n] = size(W);
0303 switch props.estmode
0304     case 'simple'
0305         for i = 1 : k
0306             S(i,1) = slevalfunctor(props.estfunctor, [], X, n, W(i, :)); 
0307         end        
0308     case 'innermul'
0309         S =  slevalfunctor(props.estfunctor, [], X, n, W, []);                      
0310 end
0311 
0312 cw = sum(W, 2);
0313 cw = cw / sum(cw);
0314 
0315 
0316 % The function to update posteriori (E-step)
0317 function [pp, W] = update_weightmap(S, cw, X, n, props)
0318 
0319 % compute conditional pdf
0320 k = length(cw);
0321 switch props.estmode
0322     case 'simple'
0323         condp = zeros(k, n);
0324         for i = 1 : k
0325             condp(i, :) = slevalfunctor(props.evalfunctor, S(i), X, n);
0326         end
0327     case 'innermul'
0328         condp = slevalfunctor(props.evalfunctor, S, X, n);
0329 end
0330 
0331 % compute posteriori
0332 switch props.condpmode
0333     case 'normal'
0334         pp = slposteriori(condp, cw);
0335     case 'log'
0336         pp = slposteriori(condp, cw, 'log');
0337 end
0338 
0339 % compute weight map
0340 W = weightmap(pp, props.weights);
0341 
0342 
0343 % The function to update models (M-step)
0344 function [S, W, cw] = update_models(S, X, W, selinds, props) 
0345 
0346 % update weights
0347 k0 = size(W, 1);
0348 [S, W, cw] = update_compweights(S, W, props);
0349 k1 = size(W, 1);
0350 if k0 ~= k1
0351     return;
0352 end
0353 
0354 % update model parameters
0355 [k, n] = size(W);
0356 switch props.estmode
0357     case 'simple'
0358         if isempty(selinds)
0359             for i = 1 : k
0360                 S(i,1) = slevalfunctor(props.estfunctor, S(i), X, n, W(i,:));
0361             end
0362         else
0363             ns = length(selinds);
0364             for ii = 1 : ns
0365                 i = selinds(ii);
0366                 S(i,1) = slevalfunctor(props.estfunctor, S(i), X, n, W(i,:));
0367             end
0368         end
0369     case 'innermul'        
0370         S = slevalfunctor(props.estfunctor, S, X, n, W, selinds);
0371 end
0372 
0373 
0374 
0375 
0376 %% Auxiliary functions
0377 
0378 % The function to select weak components
0379 
0380 function is_weak = select_weak_components(cw, thres)
0381 
0382 wthres = thres / length(cw);
0383 is_weak = (cw < wthres);
0384 
0385 % The function anneal components
0386 
0387 function [S, W, cw] = anneal_components(S, W, cw, props)
0388 
0389 is_weak = select_weak_components(cw, props.annthres);
0390 if all(is_weak)
0391     [maxw, si] = max(cw);
0392     slignorevars(maxw);
0393 elseif any(is_weak)
0394     si = find(~is_weak);
0395 else
0396     return;
0397 end
0398        
0399 switch props.estmode
0400     case 'simple'
0401         S = S(si);
0402     case 'innermul'
0403         S = feval(props.manifunc, S, 'select', si);
0404 end
0405 
0406 W = W(si, :);
0407 
0408 cw = cw(si);
0409 cw = cw / sum(cw);
0410     
0411 % The function to update weights (possibly anneal components)
0412 
0413 function [S, W, cw] = update_compweights(S, W, props)
0414 
0415 cw = sum(W, 2);
0416 cw = cw / sum(cw);
0417 
0418 if props.annthres > 0
0419     [S, W, cw] = anneal_components(S, W, cw, props);
0420 end
0421 
0422 % The function to compute weights from labels
0423 
0424 function [pp, W] = labels2weights(labels, sweights)
0425 
0426 labelset = unique(labels);
0427 [tf, inds] = ismember(labels, labelset);
0428 slignorevars(tf);
0429 clear tf;
0430 k = length(labelset);
0431 n = length(labels);
0432 
0433 pp = zeros(k, n);
0434 pp(((1:n) - 1) * k + inds) = 1;
0435 W = weightmap(pp, sweights);
0436 
0437 
0438 % The function to update posteriori to sample-component weights
0439 
0440 function W = weightmap(pp, weights)
0441 
0442 if isempty(weights)
0443     W = pp;
0444 else
0445     W = slmulvec(pp, weights, 2);
0446 end
0447 
0448 
0449 %% Simple Auxiliary Functions
0450 
0451 function checkvalid(name, val, range)
0452 
0453 if ~ismember(val, range)
0454     error('sltoolbox:invalidarg', ...
0455         'Invalid %s: %s', name, val);
0456 end
0457 
0458 
0459
Generated on Wed 20-Sep-2006 12:43:11 by m2html © 2003

Contact us at files@mathworks.com