comparison toolboxes/FullBNT-1.0.7/KPMstats/cwr_em.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function cwr  = cwr_em(X, Y, nc, varargin)
+% CWR_LEARN Fit the parameters of a cluster weighted regression model using EM
+% function cwr  = cwr_learn(X, Y, ...)
+%
+% X(:, t) is the t'th input example
+% Y(:, t) is the t'th output example
+% nc is the number of clusters
+%
+% Kevin Murphy, May 2003
+
+[max_iter, thresh, cov_typeX, cov_typeY, clamp_weights, ...
+ muX, muY, SigmaX, SigmaY, weightsY, priorC, create_init_params, ...
+cov_priorX, cov_priorY, verbose, regress, clamp_covX, clamp_covY] = process_options(...
+    varargin, 'max_iter', 10, 'thresh', 1e-2, 'cov_typeX', 'full', ...
+     'cov_typeY', 'full', 'clamp_weights', 0, ...
+     'muX', [], 'muY', [], 'SigmaX', [], 'SigmaY', [], 'weightsY', [], 'priorC', [], ...
+     'create_init_params', 1, 'cov_priorX', [], 'cov_priorY', [], 'verbose', 0, ...
+    'regress', 1, 'clamp_covX', 0, 'clamp_covY', 0);
+     
+[nx N] = size(X);
+[ny N2] = size(Y);
+if N ~= N2
+  error(sprintf('nsamples X (%d) ~= nsamples Y (%d)', N, N2));
+end
+%if N < nx 
+%  fprintf('cwr_em warning: dim X (%d) > nsamples X (%d)\n', nx, N);
+%end
+if (N < nx) & regress
+  fprintf('cwr_em warning: dim X = %d, nsamples X = %d\n', nx, N);
+end
+if (N < ny) 
+  fprintf('cwr_em warning: dim Y = %d, nsamples Y = %d\n', ny, N);
+end
+if (nc > N) 
+  error(sprintf('cwr_em: more centers (%d) than data', nc))
+end
+
+if nc==1
+  % No latent variable, so there is a closed-form solution
+  w = 1/N;
+  WYbig = Y*w;
+  WYY = WYbig * Y'; 
+  WY = sum(WYbig, 2);
+  WYTY = sum(diag(WYbig' * Y));
+  cwr.priorC = 1;
+  cwr.SigmaX = [];
+  if ~regress
+    % This is just fitting an unconditional Gaussian
+    cwr.weightsY = [];
+    [cwr.muY, cwr.SigmaY] = ...
+	mixgauss_Mstep(1, WY, WYY, WYTY, ...
+		       'cov_type', cov_typeY, 'cov_prior', cov_priorY);
+    % There is a much easier way...
+    assert(approxeq(cwr.muY, mean(Y')))
+    assert(approxeq(cwr.SigmaY, cov(Y') + 0.01*eye(ny)))
+  else
+    % This is just linear regression
+    WXbig = X*w;
+    WXX = WXbig * X';
+   WX = sum(WXbig, 2);
+    WXTX = sum(diag(WXbig' * X));
+    WXY = WXbig * Y';
+    [cwr.muY, cwr.SigmaY, cwr.weightsY] = ...
+	clg_Mstep(1, WY, WYY, WYTY, WX, WXX, WXY, ...
+		  'cov_type', cov_typeY, 'cov_prior', cov_priorY);
+  end
+  if clamp_covY, cwr.SigmaY = SigmaY; end
+  if clamp_weights,  cwr.weightsY = weightsY; end
+  return;
+end
+
+
+if create_init_params
+  [cwr.muX, cwr.SigmaX] = mixgauss_init(nc, X, cov_typeX);
+  [cwr.muY, cwr.SigmaY] = mixgauss_init(nc, Y, cov_typeY);
+  cwr.weightsY = zeros(ny, nx, nc);
+  cwr.priorC = normalize(ones(nc,1));
+else
+  cwr.muX = muX;  cwr.muY = muY; cwr.SigmaX = SigmaX; cwr.SigmaY = SigmaY;
+  cwr.weightsY = weightsY; cwr.priorC = priorC;
+end
+
+
+if clamp_covY, cwr.SigmaY = SigmaY; end
+if clamp_covX,  cwr.SigmaX = SigmaX; end
+if clamp_weights,  cwr.weightsY = weightsY; end
+
+previous_loglik = -inf;
+num_iter = 1;
+converged = 0;
+
+while (num_iter <= max_iter) & ~converged
+
+  % E step
+  
+  [likXandY, likYgivenX, post] = cwr_prob(cwr, X, Y);
+  loglik = sum(log(likXandY));
+  % extract expected sufficient statistics
+  w = sum(post,2);  % post(c,t)
+  WYY = zeros(ny, ny, nc);
+  WY = zeros(ny, nc);
+  WYTY = zeros(nc,1);
+  
+  WXX = zeros(nx, nx, nc);
+  WX = zeros(nx, nc);
+  WXTX = zeros(nc, 1);
+  WXY = zeros(nx,ny,nc);
+  %WYY = repmat(reshape(w, [1 1 nc]), [ny ny 1]) .*  repmat(Y*Y', [1 1 nc]);
+  for c=1:nc
+    weights = repmat(post(c,:), ny, 1);
+    WYbig = Y .* weights;
+    WYY(:,:,c) = WYbig * Y'; 
+    WY(:,c) = sum(WYbig, 2);
+    WYTY(c) = sum(diag(WYbig' * Y));
+
+    weights = repmat(post(c,:), nx, 1); % weights(nx, nsamples)
+    WXbig = X .* weights;
+    WXX(:,:,c) = WXbig * X';
+    WX(:,c) = sum(WXbig, 2);
+    WXTX(c) = sum(diag(WXbig' * X));
+    WXY(:,:,c) = WXbig * Y';
+  end
+
+  % M step
+  % Q -> X is called Q->Y in Mstep_clg
+  [cwr.muX, cwr.SigmaX] = mixgauss_Mstep(w, WX, WXX, WXTX, ...
+			    'cov_type', cov_typeX, 'cov_prior', cov_priorX);
+  for c=1:nc
+    assert(is_psd(cwr.SigmaX(:,:,c)))
+  end
+  
+  if clamp_weights % affects estimate of mu and Sigma
+    W = cwr.weightsY;
+  else
+    W = [];
+  end
+  [cwr.muY, cwr.SigmaY, cwr.weightsY] = ...
+      clg_Mstep(w, WY, WYY, WYTY, WX, WXX, WXY, ...
+		'cov_type', cov_typeY, 'clamped_weights', W, ...
+		'cov_prior', cov_priorY);
+  %'xs', X, 'ys', Y, 'post', post); % debug
+  %a = linspace(min(Y(2,:)), max(Y(2,:)), nc+2);
+  %cwr.muY(2,:) = a(2:end-1);
+
+  cwr.priorC = normalize(w);
+
+  for c=1:nc
+    assert(is_psd(cwr.SigmaY(:,:,c)))
+  end
+
+  if clamp_covY, cwr.SigmaY = SigmaY; end
+  if clamp_covX,  cwr.SigmaX = SigmaX; end
+  if clamp_weights,  cwr.weightsY = weightsY; end
+
+  if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
+  num_iter =  num_iter + 1;
+  converged = em_converged(loglik, previous_loglik, thresh);
+  previous_loglik = loglik;
+  
+end
+