wolffd@0: function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY)
wolffd@0: % CLG_MSTEP_SIMPLE Same as CLG_MSTEP, but doesn;t estimate Sigma,  so is slightly faster
wolffd@0: % function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY)
wolffd@0: %
wolffd@0: % See clg_Mstep for details.
wolffd@0: % Unlike clg_Mstep, there are no optional arguments, which are slow to process
wolffd@0: % if this function is inside a tight loop.
wolffd@0: 
wolffd@0: [Ysz Q] = size(Y);
wolffd@0: 
wolffd@0: if isempty(X) % no regression
wolffd@0:   %B = [];
wolffd@0:   B2 = zeros(Ysz, 1, Q);
wolffd@0:   for i=1:Q
wolffd@0:     B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q
wolffd@0:   end
wolffd@0:   [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY);
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: N = sum(w);
wolffd@0: %YY = YY + cov_prior; % regularize the scatter matrix
wolffd@0: 
wolffd@0: % Set any zero weights to one before dividing
wolffd@0: % This is valid because w(i)=0 => Y(:,i)=0, etc
wolffd@0: w = w + (w==0);
wolffd@0: 
wolffd@0: Xsz = size(X,1);
wolffd@0: % Append 1 to X to get Z
wolffd@0: ZZ = zeros(Xsz+1, Xsz+1, Q);
wolffd@0: ZY = zeros(Xsz+1, Ysz, Q);
wolffd@0: for i=1:Q
wolffd@0:   ZZ(:,:,i) = [XX(:,:,i)  X(:,i);
wolffd@0: 	       X(:,i)'    w(i)];
wolffd@0:   ZY(:,:,i) = [XY(:,:,i);
wolffd@0: 	       Y(:,i)'];
wolffd@0: end
wolffd@0: 
wolffd@0: mu = zeros(Ysz, Q);
wolffd@0: B = zeros(Ysz, Xsz, Q);
wolffd@0: for i=1:Q
wolffd@0:   % eqn 9
wolffd@0:   if rcond(ZZ(:,:,i)) < 1e-10
wolffd@0:     sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i);
wolffd@0:     %probably because there are too few cases for a high-dimensional input
wolffd@0:     ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1);
wolffd@0:   end
wolffd@0:   %A = ZY(:,:,i)' * inv(ZZ(:,:,i));
wolffd@0:   A = (ZZ(:,:,i) \ ZY(:,:,i))';
wolffd@0:   B(:,:,i) = A(:, 1:Xsz);
wolffd@0:   mu(:,i) = A(:, Xsz+1);
wolffd@0: end