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