Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/netlab3.3/gmmem.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [mix, options, errlog] = gmmem(mix, x, options)
wolffd@0 2 %GMMEM EM algorithm for Gaussian mixture model.
wolffd@0 3 %
wolffd@0 4 % Description
wolffd@0 5 % [MIX, OPTIONS, ERRLOG] = GMMEM(MIX, X, OPTIONS) uses the Expectation
wolffd@0 6 % Maximization algorithm of Dempster et al. to estimate the parameters
wolffd@0 7 % of a Gaussian mixture model defined by a data structure MIX. The
wolffd@0 8 % matrix X represents the data whose expectation is maximized, with
wolffd@0 9 % each row corresponding to a vector. The optional parameters have
wolffd@0 10 % the following interpretations.
wolffd@0 11 %
wolffd@0 12 % OPTIONS(1) is set to 1 to display error values; also logs error
wolffd@0 13 % values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then
wolffd@0 14 % only warning messages are displayed. If OPTIONS(1) is -1, then
wolffd@0 15 % nothing is displayed.
wolffd@0 16 %
wolffd@0 17 % OPTIONS(3) is a measure of the absolute precision required of the
wolffd@0 18 % error function at the solution. If the change in log likelihood
wolffd@0 19 % between two steps of the EM algorithm is less than this value, then
wolffd@0 20 % the function terminates.
wolffd@0 21 %
wolffd@0 22 % OPTIONS(5) is set to 1 if a covariance matrix is reset to its
wolffd@0 23 % original value when any of its singular values are too small (less
wolffd@0 24 % than MIN_COVAR which has the value eps). With the default value of
wolffd@0 25 % 0 no action is taken.
wolffd@0 26 %
wolffd@0 27 % OPTIONS(14) is the maximum number of iterations; default 100.
wolffd@0 28 %
wolffd@0 29 % The optional return value OPTIONS contains the final error value
wolffd@0 30 % (i.e. data log likelihood) in OPTIONS(8).
wolffd@0 31 %
wolffd@0 32 % See also
wolffd@0 33 % GMM, GMMINIT
wolffd@0 34 %
wolffd@0 35
wolffd@0 36 % Copyright (c) Ian T Nabney (1996-2001)
wolffd@0 37
wolffd@0 38 % Check that inputs are consistent
wolffd@0 39 errstring = consist(mix, 'gmm', x);
wolffd@0 40 if ~isempty(errstring)
wolffd@0 41 error(errstring);
wolffd@0 42 end
wolffd@0 43
wolffd@0 44 [ndata, xdim] = size(x);
wolffd@0 45
wolffd@0 46 % Sort out the options
wolffd@0 47 if (options(14))
wolffd@0 48 niters = options(14);
wolffd@0 49 else
wolffd@0 50 niters = 100;
wolffd@0 51 end
wolffd@0 52
wolffd@0 53 display = options(1);
wolffd@0 54 store = 0;
wolffd@0 55 if (nargout > 2)
wolffd@0 56 store = 1; % Store the error values to return them
wolffd@0 57 errlog = zeros(1, niters);
wolffd@0 58 end
wolffd@0 59 test = 0;
wolffd@0 60 if options(3) > 0.0
wolffd@0 61 test = 1; % Test log likelihood for termination
wolffd@0 62 end
wolffd@0 63
wolffd@0 64 check_covars = 0;
wolffd@0 65 if options(5) >= 1
wolffd@0 66 if display >= 0
wolffd@0 67 disp('check_covars is on');
wolffd@0 68 end
wolffd@0 69 check_covars = 1; % Ensure that covariances don't collapse
wolffd@0 70 MIN_COVAR = eps; % Minimum singular value of covariance matrix
wolffd@0 71 init_covars = mix.covars;
wolffd@0 72 end
wolffd@0 73
wolffd@0 74 % Main loop of algorithm
wolffd@0 75 for n = 1:niters
wolffd@0 76
wolffd@0 77 % Calculate posteriors based on old parameters
wolffd@0 78 [post, act] = gmmpost(mix, x);
wolffd@0 79
wolffd@0 80 % Calculate error value if needed
wolffd@0 81 if (display | store | test)
wolffd@0 82 prob = act*(mix.priors)';
wolffd@0 83 % Error value is negative log likelihood of data
wolffd@0 84 e = - sum(log(prob));
wolffd@0 85 if store
wolffd@0 86 errlog(n) = e;
wolffd@0 87 end
wolffd@0 88 if display > 0
wolffd@0 89 fprintf(1, 'Cycle %4d Error %11.6f\n', n, e);
wolffd@0 90 end
wolffd@0 91 if test
wolffd@0 92 if (n > 1 & abs(e - eold) < options(3))
wolffd@0 93 options(8) = e;
wolffd@0 94 return;
wolffd@0 95 else
wolffd@0 96 eold = e;
wolffd@0 97 end
wolffd@0 98 end
wolffd@0 99 end
wolffd@0 100
wolffd@0 101 % Adjust the new estimates for the parameters
wolffd@0 102 new_pr = sum(post, 1);
wolffd@0 103 new_c = post' * x;
wolffd@0 104
wolffd@0 105 % Now move new estimates to old parameter vectors
wolffd@0 106 mix.priors = new_pr ./ ndata;
wolffd@0 107
wolffd@0 108 mix.centres = new_c ./ (new_pr' * ones(1, mix.nin));
wolffd@0 109
wolffd@0 110 switch mix.covar_type
wolffd@0 111 case 'spherical'
wolffd@0 112 n2 = dist2(x, mix.centres);
wolffd@0 113 for j = 1:mix.ncentres
wolffd@0 114 v(j) = (post(:,j)'*n2(:,j));
wolffd@0 115 end
wolffd@0 116 mix.covars = ((v./new_pr))./mix.nin;
wolffd@0 117 if check_covars
wolffd@0 118 % Ensure that no covariance is too small
wolffd@0 119 for j = 1:mix.ncentres
wolffd@0 120 if mix.covars(j) < MIN_COVAR
wolffd@0 121 mix.covars(j) = init_covars(j);
wolffd@0 122 end
wolffd@0 123 end
wolffd@0 124 end
wolffd@0 125 case 'diag'
wolffd@0 126 for j = 1:mix.ncentres
wolffd@0 127 diffs = x - (ones(ndata, 1) * mix.centres(j,:));
wolffd@0 128 mix.covars(j,:) = sum((diffs.*diffs).*(post(:,j)*ones(1, ...
wolffd@0 129 mix.nin)), 1)./new_pr(j);
wolffd@0 130 end
wolffd@0 131 if check_covars
wolffd@0 132 % Ensure that no covariance is too small
wolffd@0 133 for j = 1:mix.ncentres
wolffd@0 134 if min(mix.covars(j,:)) < MIN_COVAR
wolffd@0 135 mix.covars(j,:) = init_covars(j,:);
wolffd@0 136 end
wolffd@0 137 end
wolffd@0 138 end
wolffd@0 139 case 'full'
wolffd@0 140 for j = 1:mix.ncentres
wolffd@0 141 diffs = x - (ones(ndata, 1) * mix.centres(j,:));
wolffd@0 142 diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
wolffd@0 143 mix.covars(:,:,j) = (diffs'*diffs)/new_pr(j);
wolffd@0 144 end
wolffd@0 145 if check_covars
wolffd@0 146 % Ensure that no covariance is too small
wolffd@0 147 for j = 1:mix.ncentres
wolffd@0 148 if min(svd(mix.covars(:,:,j))) < MIN_COVAR
wolffd@0 149 mix.covars(:,:,j) = init_covars(:,:,j);
wolffd@0 150 end
wolffd@0 151 end
wolffd@0 152 end
wolffd@0 153 case 'ppca'
wolffd@0 154 for j = 1:mix.ncentres
wolffd@0 155 diffs = x - (ones(ndata, 1) * mix.centres(j,:));
wolffd@0 156 diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin));
wolffd@0 157 [tempcovars, tempU, templambda] = ...
wolffd@0 158 ppca((diffs'*diffs)/new_pr(j), mix.ppca_dim);
wolffd@0 159 if length(templambda) ~= mix.ppca_dim
wolffd@0 160 error('Unable to extract enough components');
wolffd@0 161 else
wolffd@0 162 mix.covars(j) = tempcovars;
wolffd@0 163 mix.U(:, :, j) = tempU;
wolffd@0 164 mix.lambda(j, :) = templambda;
wolffd@0 165 end
wolffd@0 166 end
wolffd@0 167 if check_covars
wolffd@0 168 if mix.covars(j) < MIN_COVAR
wolffd@0 169 mix.covars(j) = init_covars(j);
wolffd@0 170 end
wolffd@0 171 end
wolffd@0 172 otherwise
wolffd@0 173 error(['Unknown covariance type ', mix.covar_type]);
wolffd@0 174 end
wolffd@0 175 end
wolffd@0 176
wolffd@0 177 options(8) = -sum(log(gmmprob(mix, x)));
wolffd@0 178 if (display >= 0)
wolffd@0 179 disp(maxitmess);
wolffd@0 180 end
wolffd@0 181