matthiasm@8: function [muY, SigmaY, weightsY]  = linear_regression(X, Y, varargin)
matthiasm@8: % LINEAR_REGRESSION Fit params for P(Y|X) = N(Y;  W X + mu, Sigma) 
matthiasm@8: %
matthiasm@8: % X(:, t) is the t'th input example
matthiasm@8: % Y(:, t) is the t'th output example
matthiasm@8: %
matthiasm@8: % Kevin Murphy, August 2003
matthiasm@8: %
matthiasm@8: % This is a special case of cwr_em with 1 cluster.
matthiasm@8: % You can also think of it as a front end to clg_Mstep.
matthiasm@8: 
matthiasm@8: [cov_typeY, clamp_weights,  muY, SigmaY, weightsY,...
matthiasm@8:  cov_priorY,  regress, clamp_covY] = process_options(...
matthiasm@8:     varargin, ...
matthiasm@8:      'cov_typeY', 'full', 'clamp_weights', 0, ...
matthiasm@8:      'muY', [], 'SigmaY', [], 'weightsY', [], ...
matthiasm@8:      'cov_priorY', [], 'regress', 1,  'clamp_covY', 0);
matthiasm@8:      
matthiasm@8: [nx N] = size(X);
matthiasm@8: [ny N2] = size(Y);
matthiasm@8: if N ~= N2
matthiasm@8:   error(sprintf('nsamples X (%d) ~= nsamples Y (%d)', N, N2));
matthiasm@8: end
matthiasm@8: 
matthiasm@8: w = 1/N;
matthiasm@8: WYbig = Y*w;
matthiasm@8: WYY = WYbig * Y'; 
matthiasm@8: WY = sum(WYbig, 2);
matthiasm@8: WYTY = sum(diag(WYbig' * Y));
matthiasm@8: if ~regress
matthiasm@8:   % This is just fitting an unconditional Gaussian
matthiasm@8:   weightsY = [];
matthiasm@8:   [muY, SigmaY] = ...
matthiasm@8:       mixgauss_Mstep(1, WY, WYY, WYTY, ...
matthiasm@8: 		     'cov_type', cov_typeY, 'cov_prior', cov_priorY);
matthiasm@8:   % There is a much easier way...
matthiasm@8:   assert(approxeq(muY, mean(Y')))
matthiasm@8:   assert(approxeq(SigmaY, cov(Y') + 0.01*eye(ny)))
matthiasm@8: else
matthiasm@8:   % This is just linear regression
matthiasm@8:   WXbig = X*w;
matthiasm@8:   WXX = WXbig * X';
matthiasm@8:   WX = sum(WXbig, 2);
matthiasm@8:   WXTX = sum(diag(WXbig' * X));
matthiasm@8:   WXY = WXbig * Y';
matthiasm@8:   [muY, SigmaY, weightsY] = ...
matthiasm@8:       clg_Mstep(1, WY, WYY, WYTY, WX, WXX, WXY, ...
matthiasm@8: 		'cov_type', cov_typeY, 'cov_prior', cov_priorY);
matthiasm@8: end
matthiasm@8: if clamp_covY, SigmaY = SigmaY; end
matthiasm@8: if clamp_weights,  weightsY = weightsY; end
matthiasm@8: 
matthiasm@8: if nx==1 & ny==1 & regress
matthiasm@8:   P = polyfit(X,Y); % Y = P(1) X^1 + P(2) X^0 = ax + b
matthiasm@8:   assert(approxeq(muY, P(2)))
matthiasm@8:   assert(approxeq(weightsY, P(1)))
matthiasm@8: end
matthiasm@8: 
matthiasm@8: %%%%%%%% Test
matthiasm@8: if 0 
matthiasm@8:   c1 = randn(2,100);   c2 = randn(2,100);
matthiasm@8:   y = c2(1,:); X = [ones(size(c1,2),1) c1'];
matthiasm@8:   b = regress(y(:), X); % stats toolbox
matthiasm@8:   [m,s,w] = linear_regression(c1, y);
matthiasm@8:   assert(approxeq(b(1),m))
matthiasm@8:   assert(approxeq(b(2), w(1)))
matthiasm@8:   assert(approxeq(b(3), w(2)))
matthiasm@8: end