clg_mstep_simple.m~

matlab Netlab库

function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY)
% CLG_MSTEP_SIMPLE Same as CLG_MSTEP, but doesn;t estimate Sigma,  so is slightly faster
% function [mu, B] = clg_Mstep(w, Y, YY, YTY, X, XX, XY)
%
% We fit P(Y|X,Q=i) = N(Y; B_i X + mu_i, Sigma_i) 
% and w(i,t) = p(M(t)=i|y(t)) = posterior responsibility
% See www.ai.mit.edu/~murphyk/Papers/learncg.pdf.
%
%
% INPUTS:
% w(i) = sum_t w(i,t) = responsibilities for each mixture component
%  If there is only one mixture component (i.e., Q does not exist),
%  then w(i) = N = nsamples,  and 
%  all references to i can be replaced by 1.
% Y(:,i) = sum_t w(i,t) y(:,t) = weighted observations
% YY(:,:,i) = sum_t w(i,t) y(:,t) y(:,t)' = weighted outer product
% YTY(i) = sum_t w(i,t) y(:,t)' y(:,t) = weighted inner product
%   You only need to pass in YTY if Sigma is to be estimated as spherical.
%
% In the regression context, we must also pass in the following
% X(:,i) = sum_t w(i,t) x(:,t) = weighted inputs
% XX(:,:,i) = sum_t w(i,t) x(:,t) x(:,t)' = weighted outer product
% XY(i) = sum_t w(i,t) x(:,t) y(:,t)' = weighted outer product
%
% Unlike clg_Mstep, there are no optional arguments, which are slow to process
% if this function is inside a tight loop

[Ysz Q] = size(Y);

if isempty(X) % no regression
  %B = [];
  B2 = zeros(Ysz, 1, Q);
  for i=1:Q
    B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q
  end
  [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin{:});
  return;
end

N = sum(w);
%YY = YY + cov_prior; % regularize the scatter matrix

% Set any zero weights to one before dividing
% This is valid because w(i)=0 => Y(:,i)=0, etc
w = w + (w==0);

Xsz = size(X,1);
% Append 1 to X to get Z
ZZ = zeros(Xsz+1, Xsz+1, Q);
ZY = zeros(Xsz+1, Ysz, Q);
for i=1:Q
  ZZ(:,:,i) = [XX(:,:,i)  X(:,i);
	       X(:,i)'    w(i)];
  ZY(:,:,i) = [XY(:,:,i);
	       Y(:,i)'];
end

mu = zeros(Ysz, Q);
B = zeros(Ysz, Xsz, Q);
for i=1:Q
  % eqn 9
  if rcond(ZZ(:,:,i)) < 1e-10
    sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i);
    %probably because there are too few cases for a high-dimensional input
    ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1);
  end
  %A = ZY(:,:,i)' * inv(ZZ(:,:,i));
  A = (ZZ(:,:,i) \ ZY(:,:,i))';
  B(:,:,i) = A(:, 1:Xsz);
  mu(:,i) = A(:, Xsz+1);
end