V_GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,g,f]=(x,c,l,m0,v0,w0,wx)
Usage:
(1) [m,v,w]=v_gaussmix(x,[],[],k); % create GMM with k mixtures and diagonal covariances
(2) [m,v,w]=gaussmix(x,[],[],k,'v'); % create GMM with k mixtures and full covariances
Inputs: n data values, k mixtures, p parameters, l loops
X(n,p) Input data vectors, one per row.
C(1) Minimum variance of normalized data (Use [] to take default value of n^-2)
L The integer portion of l gives a maximum loop count. The fractional portion gives
an optional stopping threshold. Iteration will cease if the increase in
log likelihood density per data point is less than this value. Thus l=10.001 will
stop after 10 iterations or when the increase in log likelihood falls below
0.001.
As a special case, if L=0, then the first three outputs are omitted.
Use [] to take default value of 100.0001
M0 Number of mixtures required (or initial mixture means, M0(k,p), - see below)
V0 Initialization mode:
'f' Initialize with K randomly selected data points [default]
'p' Initialize with centroids and variances of random partitions
'k' k-means algorithm ('kf' and 'kp' determine initialization)
'h' k-harmonic means algorithm ('hf' and 'hp' determine initialization) [default]
's' use unscaled data during initialization phase instead of scaling it first
'm' M0(k,p) contains the initial centres
'v' full covariance matrices
Mode 'hf' [the default] generally gives the best results but 'f' is faster and often OK
W0(n,1) Data point weights (need not sum to unity)
Alternatively, initial values for M0, V0 and W0 can be given explicitly:
M0(k,p) Initial mixture means, one row per mixture.
V0(k,p) Initial mixture variances, one row per mixture.
or V0(p,p,k) one full-covariance matrix per mixture
W0(k,1) Initial mixture weights, one per mixture. The weights must sum to unity.
WX(n,1) Data point weights (need not sum to unity)
Outputs: (Note that M, V and W are omitted if L==0)
M(k,p) Mixture means, one row per mixture. (omitted if L==0)
V(k,p) Mixture variances, one row per mixture. (omitted if L==0)
or V(p,p,k) if full covariance matrices in use (i.e. either 'v' option or V0(p,p,k) specified)
W(k,1) Mixture weights, one per mixture. The weights will sum to unity. (omitted if L==0)
G Average log probability of the input data points.
F Fisher's Discriminant measures how well the data divides into classes.
It is the ratio of the between-mixture variance to the average mixture variance: a
high value means the classes (mixtures) are well separated.
PP(n,1) Log probability of each data point
GG(l+1,1) Average log probabilities at the beginning of each iteration and at the end of the last
The fitting procedure uses one of several initialization methods to create an initial guess
for the mixture centres and then uses the EM (expectation-maximization) algorithm to refine
the guess. Although the EM algorithm is deterministic, the initialization procedures use
random numbers and so the routine will not give identical answers if you call it multiple
times with the same input data. See v_randvec() for generating GMM data vectors.0001 function [m,v,w,g,f,pp,gg]=v_gaussmix(x,c,l,m0,v0,w0,wx) 0002 %V_GAUSSMIX fits a gaussian mixture pdf to a set of data observations [m,v,w,g,f]=(x,c,l,m0,v0,w0,wx) 0003 % 0004 % Usage: 0005 % (1) [m,v,w]=v_gaussmix(x,[],[],k); % create GMM with k mixtures and diagonal covariances 0006 % (2) [m,v,w]=gaussmix(x,[],[],k,'v'); % create GMM with k mixtures and full covariances 0007 % 0008 % Inputs: n data values, k mixtures, p parameters, l loops 0009 % 0010 % X(n,p) Input data vectors, one per row. 0011 % C(1) Minimum variance of normalized data (Use [] to take default value of n^-2) 0012 % L The integer portion of l gives a maximum loop count. The fractional portion gives 0013 % an optional stopping threshold. Iteration will cease if the increase in 0014 % log likelihood density per data point is less than this value. Thus l=10.001 will 0015 % stop after 10 iterations or when the increase in log likelihood falls below 0016 % 0.001. 0017 % As a special case, if L=0, then the first three outputs are omitted. 0018 % Use [] to take default value of 100.0001 0019 % M0 Number of mixtures required (or initial mixture means, M0(k,p), - see below) 0020 % V0 Initialization mode: 0021 % 'f' Initialize with K randomly selected data points [default] 0022 % 'p' Initialize with centroids and variances of random partitions 0023 % 'k' k-means algorithm ('kf' and 'kp' determine initialization) 0024 % 'h' k-harmonic means algorithm ('hf' and 'hp' determine initialization) [default] 0025 % 's' use unscaled data during initialization phase instead of scaling it first 0026 % 'm' M0(k,p) contains the initial centres 0027 % 'v' full covariance matrices 0028 % Mode 'hf' [the default] generally gives the best results but 'f' is faster and often OK 0029 % W0(n,1) Data point weights (need not sum to unity) 0030 % 0031 % Alternatively, initial values for M0, V0 and W0 can be given explicitly: 0032 % 0033 % M0(k,p) Initial mixture means, one row per mixture. 0034 % V0(k,p) Initial mixture variances, one row per mixture. 0035 % or V0(p,p,k) one full-covariance matrix per mixture 0036 % W0(k,1) Initial mixture weights, one per mixture. The weights must sum to unity. 0037 % WX(n,1) Data point weights (need not sum to unity) 0038 % 0039 % Outputs: (Note that M, V and W are omitted if L==0) 0040 % 0041 % M(k,p) Mixture means, one row per mixture. (omitted if L==0) 0042 % V(k,p) Mixture variances, one row per mixture. (omitted if L==0) 0043 % or V(p,p,k) if full covariance matrices in use (i.e. either 'v' option or V0(p,p,k) specified) 0044 % W(k,1) Mixture weights, one per mixture. The weights will sum to unity. (omitted if L==0) 0045 % G Average log probability of the input data points. 0046 % F Fisher's Discriminant measures how well the data divides into classes. 0047 % It is the ratio of the between-mixture variance to the average mixture variance: a 0048 % high value means the classes (mixtures) are well separated. 0049 % PP(n,1) Log probability of each data point 0050 % GG(l+1,1) Average log probabilities at the beginning of each iteration and at the end of the last 0051 % 0052 % The fitting procedure uses one of several initialization methods to create an initial guess 0053 % for the mixture centres and then uses the EM (expectation-maximization) algorithm to refine 0054 % the guess. Although the EM algorithm is deterministic, the initialization procedures use 0055 % random numbers and so the routine will not give identical answers if you call it multiple 0056 % times with the same input data. See v_randvec() for generating GMM data vectors. 0057 0058 % Bugs/Suggestions 0059 % (1) Allow processing in chunks by outputting/reinputting an array of sufficient statistics 0060 % (2) Other initialization options: 0061 % 'l' LBG algorithm 0062 % 'm' Move-means (dog-rabbit) algorithm 0063 % (3) Allow updating of weights-only, not means/variances 0064 % (4) Allow freezing of means and/or variances 0065 % (5) Allow initialization of means only 0066 0067 % Copyright (C) Mike Brookes 2000-2009 0068 % Version: $Id: v_gaussmix.m 10865 2018-09-21 17:22:45Z dmb $ 0069 % 0070 % VOICEBOX is a MATLAB toolbox for speech processing. 0071 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0072 % 0073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0074 % This program is free software; you can redistribute it and/or modify 0075 % it under the terms of the GNU Lesser General Public License as published by 0076 % the Free Software Foundation; either version 3 of the License, or 0077 % (at your option) any later version. 0078 % 0079 % This program is distributed in the hope that it will be useful, 0080 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0081 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0082 % GNU Lesser General Public License for more details. 0083 % 0084 % You can obtain a copy of the GNU Lesser General Public License from 0085 % https://www.gnu.org/licenses/ . 0086 % See files gpl-3.0.txt and lgpl-3.0.txt included in this distribution. 0087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0088 [n,p]=size(x); % n = number of training values, p = dimension of data vector 0089 wn=ones(n,1); 0090 memsize=v_voicebox('memsize'); % set memory size to use 0091 if isempty(c) % no minimum variance is specified 0092 c=1/n^2; 0093 else 0094 c=c(1); % just to prevent legacy code failing 0095 end 0096 fulliv=0; % initial variance is not full 0097 if isempty(l) 0098 l=100+1e-4; % max loop count + stopping threshold 0099 end 0100 if nargin<5 || isempty(v0) || ischar(v0) % no initial values specified for m0, v0, w0 0101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0102 % No initial values given, so we must use k-means or equivalent % 0103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0104 if nargin<6 0105 if nargin<5 || isempty(v0) 0106 v0='hf'; % default initialization mode: hf 0107 end 0108 wx=wn; % no data point weights 0109 else 0110 wx=w0(:); % data point weights 0111 end 0112 wx=wx(:)/sum(wx); % normalize and force to be a column vector 0113 if any(v0=='m') 0114 k=size(m0,1); 0115 else 0116 k=m0; 0117 end 0118 fv=any(v0=='v'); % full covariance matrices requested 0119 mx0=wx'*x; % calculate mean of input data in each dimension 0120 vx0=wx'*x.^2-mx0.^2; % calculate variance of input data in each dimension 0121 sx0=sqrt(vx0)+(vx0==0); % scale factor: std of input data in each dimension (or 1 if std==0) 0122 if n<=k % each data point can have its own mixture 0123 xs=(x-mx0(wn,:))./sx0(wn,:); % scale the data 0124 m=xs(mod((1:k)-1,n)+1,:); % just include all points several times 0125 v=zeros(k,p); % will be set to floor later 0126 w=zeros(k,1); 0127 w(1:n)=1/n; 0128 if l>0 0129 l=0.1; % no point in iterating 0130 end 0131 else % more points than mixtures 0132 if any(v0=='s') 0133 xs=x; % do not scale data during initialization 0134 else 0135 xs=(x-mx0(wn,:))./sx0(wn,:); % else scale now 0136 if any(v0=='m') 0137 m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:); % scale specified means as well 0138 end 0139 end 0140 w=repmat(1/k,k,1); % all mixtures equally likely 0141 if any(v0=='k') % k-means initialization 0142 if any(v0=='m') 0143 [m,e,j]=v_kmeans(xs,k,m); 0144 elseif any(v0=='p') 0145 [m,e,j]=v_kmeans(xs,k,'p'); 0146 else 0147 [m,e,j]=v_kmeans(xs,k,'f'); 0148 end 0149 elseif any(v0=='h') % k-harmonic means initialization 0150 if any(v0=='m') 0151 [m,e,j]=v_kmeanhar(xs,k,[],4,m); 0152 else 0153 if any(v0=='p') 0154 [m,e,j]=v_kmeanhar(xs,k,[],4,'p'); 0155 else 0156 [m,e,j]=v_kmeanhar(xs,k,[],4,'f'); 0157 end 0158 end 0159 elseif any(v0=='p') % Initialize using a random partition 0160 j=ceil(rand(n,1)*k); % allocate to random clusters 0161 j(v_rnsubset(k,n))=1:k; % but force at least one point per cluster 0162 for i=1:k 0163 m(i,:)=mean(xs(j==i,:),1); 0164 end 0165 else 0166 if any(v0=='m') 0167 m=m0; % use specified centres 0168 else 0169 m=xs(v_rnsubset(k,n),:); % Forgy initialization: sample k centres without replacement [default] 0170 end 0171 [e,j]=v_kmeans(xs,k,m,0); % find out the cluster allocation 0172 end 0173 if any(v0=='s') 0174 xs=(x-mx0(wn,:))./sx0(wn,:); % scale data now if not done previously 0175 end 0176 v=zeros(k,p); % initialise diagonal covariances 0177 w=zeros(k,1); % initialise mixture weights 0178 for i=1:k 0179 ni=sum(j==i); % number assigned to this centre 0180 w(i)=(ni+1)/(n+k); % weight of this mixture 0181 if ni % if any data points are assined to this mixture 0182 v(i,:)=sum((xs(j==i,:)-repmat(m(i,:),ni,1)).^2,1)/ni; % initialize mixture variance 0183 else % else if no points assigned to this mixture 0184 v(i,:)=zeros(1,p); % set mixture variance to zero (will be set to the global floor later) 0185 end 0186 end 0187 end 0188 else 0189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0190 % use initial values given as input parameters % 0191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0192 if nargin<7 0193 wx=wn; % no data point weights specified so use uniform weights 0194 end 0195 wx=wx(:)/sum(wx); % normalize weights and force a column vector 0196 mx0=wx'*x; % calculate mean of input data in each dimension 0197 vx0=wx'*x.^2-mx0.^2; % calculate variance of input data in each dimension 0198 sx0=sqrt(vx0); 0199 sx0(sx0==0)=1; % do not divide by zero when scaling 0200 [k,p]=size(m0); 0201 xs=(x-mx0(wn,:))./sx0(wn,:); % scale the data 0202 m=(m0-mx0(ones(k,1),:))./sx0(ones(k,1),:); % and the means 0203 v=v0; 0204 w=w0; 0205 fv=ndims(v)>2 || size(v,1)>k; % full covariance matrix is supplied: v(p,p,k) 0206 if fv % if full covariance matrix is supplied 0207 mk=eye(p)==0; % create mask for off-diagonal elements 0208 fulliv=any(v(repmat(mk,[1 1 k]))~=0); % check if any are non-zero 0209 if ~fulliv % supplied full covariance matrices are actually diagonal 0210 v=reshape(v(repmat(~mk,[1 1 k])),p,k)'./repmat(sx0.^2,k,1); % just pick out and scale the diagonal elements for now 0211 else 0212 v=v./repmat(sx0'*sx0,[1 1 k]); % scale the full covariance matrix for each mixture 0213 end 0214 else % diagonal covariance matrix is supplied: v(k,p) 0215 v = v ./ repmat(sx0.^2,k,1); % scale the diagonal covariance matrix for each mixture 0216 end 0217 end 0218 if length(wx)~=n 0219 error('%d datapoints but %d weights',n,length(wx)); 0220 end 0221 lsx=sum(log(sx0)); 0222 xsw=xs.*repmat(wx,1,p); % weighted data points 0223 if ~fulliv % initializing with diagonal covariance if v0 is either unspecified or diagonal 0224 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0225 % Diagonal Covariance matrices % 0226 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0227 v=max(v,c); % apply the lower bound 0228 xs2=xs.^2.*repmat(wx,1,p); % square and weight the data for variance calculations 0229 0230 % If data size is large then do calculations in chunks 0231 0232 nb=min(n,max(1,floor(memsize/(8*p*k)))); % chunk size for testing data points 0233 nl=ceil(n/nb); % number of chunks 0234 jx0=n-(nl-1)*nb; % size of first chunk 0235 0236 im=repmat(1:k,1,nb); im=im(:); 0237 th=(l-floor(l))*n; 0238 sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values 0239 lp=floor(l)+sd; % extra loop needed to calculate final G value 0240 0241 lpx=zeros(1,n); % log probability of each data point 0242 wk=ones(k,1); 0243 wp=ones(1,p); 0244 wnb=ones(1,nb); 0245 wnj=ones(1,jx0); 0246 0247 % EM loop 0248 0249 g=0; % dummy initial value for comparison 0250 gg=zeros(lp+1,1); 0251 ss=sd; % initialize stopping count (0 or 1) 0252 for j=1:lp 0253 g1=g; % save previous log likelihood (2*pi factor omitted) 0254 m1=m; % save previous means, variances and weights 0255 v1=v; 0256 w1=w; 0257 vi=-0.5*v.^(-1); % data-independent scale factor in exponent 0258 lvm=log(w)-0.5*sum(log(v),2); % log of external scale factor (excluding -0.5*p*log(2pi) term) 0259 0260 % first do partial chunk (of length jx0) 0261 0262 jx=jx0; 0263 ii=1:jx; % indices of data points in this chunk 0264 kk=repmat(ii,k,1); % kk(k,jx): data point index 0265 km=repmat(1:k,1,jx); % km(k,jx): mixture index 0266 py=reshape(sum((xs(kk(:),:)-m(km(:),:)).^2.*vi(km(:),:),2),k,jx)+lvm(:,wnj); % py(k,jx) pdf of each point with each mixture 0267 mx=max(py,[],1); % mx(1,jx) find normalizing factor for each data point to prevent underflow when using exp() 0268 px=exp(py-mx(wk,:)); % find normalized probability of each mixture for each datapoint 0269 ps=sum(px,1); % total normalized likelihood of each data point 0270 px=px./ps(wk,:); % relative mixture probabilities for each data point (columns sum to 1) 0271 lpx(ii)=log(ps)+mx; 0272 pk=px*wx(ii); % pk(k,1) effective fraction of data points for each mixture (could be zero due to underflow) 0273 sx=px*xsw(ii,:); 0274 sx2=px*xs2(ii,:); 0275 for il=2:nl % process the data points in chunks 0276 ix=jx+1; 0277 jx=jx+nb; % increment upper limit 0278 ii=ix:jx; % indices of data points in this chunk 0279 kk=repmat(ii,k,1); 0280 py=reshape(sum((xs(kk(:),:)-m(im,:)).^2.*vi(im,:),2),k,nb)+lvm(:,wnb); 0281 mx=max(py,[],1); % find normalizing factor for each data point to prevent underflow when using exp() 0282 px=exp(py-mx(wk,:)); % find normalized probability of each mixture for each datapoint 0283 ps=sum(px,1); % total normalized likelihood of each data point 0284 px=px./ps(wk,:); % relative mixture probabilities for each data point (columns sum to 1) 0285 lpx(ii)=log(ps)+mx; 0286 pk=pk+px*wx(ii); % pk(k,1) effective fraction of data points for each mixture (could be zero due to underflow) 0287 sx=sx+px*xsw(ii,:); 0288 sx2=sx2+px*xs2(ii,:); 0289 end 0290 g=lpx*wx; % total log probability summed over all data points 0291 gg(j)=g; % save log prob at each iteration 0292 w=pk; % normalize to get the weights 0293 if pk % if all elements of pk are non-zero 0294 m=sx./pk(:,wp); % calculate mixture means 0295 v=sx2./pk(:,wp); % and raw second moments 0296 else 0297 wm=pk==0; % mask indicating mixtures with zero weights 0298 nz=sum(wm); % number of zero-weight mixtures 0299 [vv,mk]=sort(lpx); % find the lowest probability data points 0300 m=zeros(k,p); % initialize means and variances to zero (variances are floored later) 0301 v=m; 0302 m(wm,:)=xs(mk(1:nz),:); % set zero-weight mixture means to worst-fitted data points 0303 w(wm)=1/n; % set these weights non-zero 0304 w=w*n/(n+nz); % normalize so the weights sum to unity 0305 wm=~wm; % mask for non-zero weights 0306 m(wm,:)=sx(wm,:)./pk(wm,wp); % recalculate means and raw second moments for mixtures with a non-zero weight 0307 v(wm,:)=sx2(wm,:)./pk(wm,wp); 0308 end 0309 v=max(v-m.^2,c); % convert raw second moments to variances and apply floor 0310 if g-g1<=th && j>1 0311 if ~ss, break; end % stop 0312 ss=ss-1; % stop next time 0313 end 0314 end 0315 if sd && ~fv % we need to calculate the final probabilities 0316 pp=lpx'-0.5*p*log(2*pi)-lsx; % log of total probability of each data point 0317 gg=gg(1:j)-0.5*p*log(2*pi)-lsx; % average log prob at each iteration 0318 g=gg(end); 0319 m=m1; % back up to previous iteration 0320 v=v1; 0321 w=w1; 0322 mm=sum(m,1)/k; 0323 f=(m(:)'*m(:)-k*mm(:)'*mm(:))/sum(v(:)); 0324 end 0325 if ~fv 0326 m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:); % unscale means 0327 v=v.*repmat(sx0.^2,k,1); % and variances 0328 else 0329 v1=v; 0330 v=zeros(p,p,k); 0331 mk=eye(p)==1; % mask for diagonal elements 0332 v(repmat(mk,[1 1 k]))=v1'; % set from v1 0333 end 0334 end 0335 if fv % check if full covariance matrices were requested 0336 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0337 % Full Covariance matrices % 0338 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0339 pl=p*(p+1)/2; % number of non-trivial elements in lower triangular matrix 0340 lix=1:p^2; 0341 cix=repmat(1:p,p,1); 0342 rix=cix'; 0343 lix(cix>rix)=[]; % index of lower triangular elements 0344 cix=cix(lix); % index of lower triangular columns 0345 rix=rix(lix); % index of lower triangular rows 0346 % dix=find(rix==cix); % index of lower triangular diagonal elements [unused] 0347 lixi=zeros(p,p); 0348 lixi(lix)=1:pl; 0349 lixi=lixi'; 0350 lixi(lix)=1:pl; % reverse index to build full matrices 0351 v=reshape(v,p^2,k); 0352 v=v(lix,:)'; % lower triangular in rows (k,p*(p+1)/2) 0353 0354 % If data size is large then do calculations in chunks 0355 0356 nb=min(n,max(1,floor(memsize/(24*p*k)))); % chunk size for testing data points 0357 nl=ceil(n/nb); % number of chunks 0358 jx0=n-(nl-1)*nb; % size of first chunk 0359 % 0360 th=(l-floor(l))*n; 0361 sd=(nargout > 3*(l~=0)); % = 1 if we are outputting log likelihood values 0362 lp=floor(l)+sd; % extra loop needed to calculate final G value 0363 % 0364 lpx=zeros(1,n); % log probability of each data point 0365 wk=ones(k,1); 0366 wp=ones(1,p); 0367 wpl=ones(1,pl); % 1 index for lower triangular matrix 0368 wnb=ones(1,nb); 0369 wnj=ones(1,jx0); 0370 0371 % EM loop 0372 0373 g=0; % dummy initial value for comparison 0374 gg=zeros(lp+1,1); 0375 ss=sd; % initialize stopping count (0 or 1) 0376 vi=zeros(p*k,p); % stack of k inverse cov matrices each size p*p 0377 vim=zeros(p*k,1); % stack of k vectors of the form inv(v)*m 0378 mtk=vim; % stack of k vectors of the form m 0379 lvm=zeros(k,1); 0380 wpk=repmat((1:p)',k,1); 0381 for j=1:lp 0382 g1=g; % save previous log likelihood (2*pi factor omitted) 0383 m1=m; % save previous means, variances and weights 0384 v1=v; 0385 w1=w; 0386 for ik=1:k 0387 % these lines added for debugging only 0388 % vk=reshape(v(k,lixi),p,p); 0389 % condk(ik)=cond(vk); 0390 %%%%%%%%%%%%%%%%%%%% 0391 [uvk,dvk]=eig(reshape(v(ik,lixi),p,p)); % convert lower triangular to full and find eigenvalues 0392 dvk=max(diag(dvk),c); % apply variance floor to eigenvalues 0393 vik=-0.5*uvk*diag(dvk.^(-1))*uvk'; % calculate inverse 0394 vi((ik-1)*p+(1:p),:)=vik; % vi contains all mixture inverses stacked on top of each other 0395 vim((ik-1)*p+(1:p))=vik*m(ik,:)'; % vim contains vi*m for all mixtures stacked on top of each other 0396 mtk((ik-1)*p+(1:p))=m(ik,:)'; % mtk contains all mixture means stacked on top of each other 0397 lvm(ik)=log(w(ik))-0.5*sum(log(dvk)); % vm contains the weighted sqrt of det(vi) for each mixture 0398 end 0399 % 0400 % % first do partial chunk 0401 % 0402 jx=jx0; 0403 ii=1:jx; 0404 xii=xs(ii,:).'; 0405 py=reshape(sum(reshape((vi*xii-vim(:,wnj)).*(xii(wpk,:)-mtk(:,wnj)),p,jx*k),1),k,jx)+lvm(:,wnj); 0406 mx=max(py,[],1); % find normalizing factor for each data point to prevent underflow when using exp() 0407 px=exp(py-mx(wk,:)); % find normalized probability of each mixture for each datapoint 0408 ps=sum(px,1); % total normalized likelihood of each data point 0409 px=px./ps(wk,:); % relative mixture probabilities for each data point (columns sum to 1) 0410 lpx(ii)=log(ps)+mx; 0411 pk=px*wx(ii); % effective fraction of data points for each mixture (could be zero due to underflow) 0412 sx=px*xsw(ii,:); 0413 sx2=px*(xsw(ii,rix).*xs(ii,cix)); % accumulator for variance calculation (lower tri cov matrix as a row) 0414 for il=2:nl 0415 ix=jx+1; 0416 jx=jx+nb; % increment upper limit 0417 ii=ix:jx; 0418 xii=xs(ii,:).'; 0419 py=reshape(sum(reshape((vi*xii-vim(:,wnb)).*(xii(wpk,:)-mtk(:,wnb)),p,nb*k),1),k,nb)+lvm(:,wnb); 0420 mx=max(py,[],1); % find normalizing factor for each data point to prevent underflow when using exp() 0421 px=exp(py-mx(wk,:)); % find normalized probability of each mixture for each datapoint 0422 ps=sum(px,1); % total normalized likelihood of each data point 0423 px=px./ps(wk,:); % relative mixture probabilities for each data point (columns sum to 1) 0424 lpx(ii)=log(ps)+mx; 0425 pk=pk+px*wx(ii); % effective fraction of data points for each mixture (could be zero due to underflow) 0426 sx=sx+px*xsw(ii,:); % accumulator for mean calculation 0427 sx2=sx2+px*(xsw(ii,rix).*xs(ii,cix)); % accumulator for variance calculation 0428 end 0429 g=lpx*wx; % total log probability summed over all data points 0430 gg(j)=g; % save convergence history 0431 w=pk; % w(k,1) normalize to get the column of weights 0432 if pk % if all elements of pk are non-zero 0433 m=sx./pk(:,wp); % find mean and mean square 0434 v=sx2./pk(:,wpl); 0435 else 0436 wm=pk==0; % mask indicating mixtures with zero weights 0437 nz=sum(wm); % number of zero-weight mixtures 0438 [vv,mk]=sort(lpx); % find the lowest probability data points 0439 m=zeros(k,p); % initialize means and variances to zero (variances are floored later) 0440 v=zeros(k,pl); 0441 m(wm,:)=xs(mk(1:nz),:); % set zero-weight mixture means to worst-fitted data points 0442 w(wm)=1/n; % set these weights non-zero 0443 w=w*n/(n+nz); % normalize so the weights sum to unity 0444 wm=~wm; % mask for non-zero weights 0445 m(wm,:)=sx(wm,:)./pk(wm,wp); % recalculate means and variances for mixtures with a non-zero weight 0446 v(wm,:)=sx2(wm,:)./pk(wm,wpl); 0447 end 0448 v=v-m(:,cix).*m(:,rix); % subtract off mean squared 0449 if g-g1<=th && j>1 0450 if ~ss, break; end % stop 0451 ss=ss-1; % stop next time 0452 end 0453 end 0454 if sd % we need to calculate the final probabilities 0455 pp=lpx'-0.5*p*log(2*pi)-lsx; % log of total probability of each data point 0456 gg=gg(1:j)-0.5*p*log(2*pi)-lsx; % average log prob at each iteration 0457 g=gg(end); 0458 % gg' % *** DEBUG ONLY *** 0459 m=m1; % back up to previous iteration 0460 v=zeros(p,p,k); % reserve spave for k full covariance matrices 0461 trv=0; % sum of variance matrix traces 0462 for ik=1:k % loop for each mixture to apply variance floor 0463 [uvk,dvk]=eig(reshape(v1(ik,lixi),p,p)); % convert lower triangular to full and find eigenvectors 0464 dvk=max(diag(dvk),c); % apply variance floor to eigenvalues 0465 v(:,:,ik)=uvk*diag(dvk)*uvk'; % reconstitute full matrix 0466 trv=trv+sum(dvk); % add trace to the sum 0467 end 0468 w=w1; 0469 mm=sum(m,1)/k; 0470 f=(m(:)'*m(:)-k*mm(:)'*mm(:))/trv; 0471 else 0472 v1=v; % lower triangular form 0473 v=zeros(p,p,k); % reserve space for k full covariance matrices 0474 for ik=1:k % loop for each mixture to apply variance floor 0475 [uvk,dvk]=eig(reshape(v1(ik,lixi),p,p)); % convert lower triangular to full and find eigenvectors 0476 dvk=max(diag(dvk),c); % apply variance floor 0477 v(:,:,ik)=uvk*diag(dvk)*uvk'; % reconstitute full matrix 0478 end 0479 end 0480 m=m.*sx0(ones(k,1),:)+mx0(ones(k,1),:); % unscale means 0481 v=v.*repmat(sx0'*sx0,[1 1 k]); 0482 end 0483 if l==0 % suppress the first three output arguments if l==0 0484 m=g; 0485 v=f; 0486 w=pp; 0487 end 0488