|
wolffd@0
|
1 function [muY, SigmaY, weightsY] = linear_regression(X, Y, varargin)
|
|
wolffd@0
|
2 % LINEAR_REGRESSION Fit params for P(Y|X) = N(Y; W X + mu, Sigma)
|
|
wolffd@0
|
3 %
|
|
wolffd@0
|
4 % X(:, t) is the t'th input example
|
|
wolffd@0
|
5 % Y(:, t) is the t'th output example
|
|
wolffd@0
|
6 %
|
|
wolffd@0
|
7 % Kevin Murphy, August 2003
|
|
wolffd@0
|
8 %
|
|
wolffd@0
|
9 % This is a special case of cwr_em with 1 cluster.
|
|
wolffd@0
|
10 % You can also think of it as a front end to clg_Mstep.
|
|
wolffd@0
|
11
|
|
wolffd@0
|
12 [cov_typeY, clamp_weights, muY, SigmaY, weightsY,...
|
|
wolffd@0
|
13 cov_priorY, regress, clamp_covY] = process_options(...
|
|
wolffd@0
|
14 varargin, ...
|
|
wolffd@0
|
15 'cov_typeY', 'full', 'clamp_weights', 0, ...
|
|
wolffd@0
|
16 'muY', [], 'SigmaY', [], 'weightsY', [], ...
|
|
wolffd@0
|
17 'cov_priorY', [], 'regress', 1, 'clamp_covY', 0);
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 [nx N] = size(X);
|
|
wolffd@0
|
20 [ny N2] = size(Y);
|
|
wolffd@0
|
21 if N ~= N2
|
|
wolffd@0
|
22 error(sprintf('nsamples X (%d) ~= nsamples Y (%d)', N, N2));
|
|
wolffd@0
|
23 end
|
|
wolffd@0
|
24
|
|
wolffd@0
|
25 w = 1/N;
|
|
wolffd@0
|
26 WYbig = Y*w;
|
|
wolffd@0
|
27 WYY = WYbig * Y';
|
|
wolffd@0
|
28 WY = sum(WYbig, 2);
|
|
wolffd@0
|
29 WYTY = sum(diag(WYbig' * Y));
|
|
wolffd@0
|
30 if ~regress
|
|
wolffd@0
|
31 % This is just fitting an unconditional Gaussian
|
|
wolffd@0
|
32 weightsY = [];
|
|
wolffd@0
|
33 [muY, SigmaY] = ...
|
|
wolffd@0
|
34 mixgauss_Mstep(1, WY, WYY, WYTY, ...
|
|
wolffd@0
|
35 'cov_type', cov_typeY, 'cov_prior', cov_priorY);
|
|
wolffd@0
|
36 % There is a much easier way...
|
|
wolffd@0
|
37 assert(approxeq(muY, mean(Y')))
|
|
wolffd@0
|
38 assert(approxeq(SigmaY, cov(Y') + 0.01*eye(ny)))
|
|
wolffd@0
|
39 else
|
|
wolffd@0
|
40 % This is just linear regression
|
|
wolffd@0
|
41 WXbig = X*w;
|
|
wolffd@0
|
42 WXX = WXbig * X';
|
|
wolffd@0
|
43 WX = sum(WXbig, 2);
|
|
wolffd@0
|
44 WXTX = sum(diag(WXbig' * X));
|
|
wolffd@0
|
45 WXY = WXbig * Y';
|
|
wolffd@0
|
46 [muY, SigmaY, weightsY] = ...
|
|
wolffd@0
|
47 clg_Mstep(1, WY, WYY, WYTY, WX, WXX, WXY, ...
|
|
wolffd@0
|
48 'cov_type', cov_typeY, 'cov_prior', cov_priorY);
|
|
wolffd@0
|
49 end
|
|
wolffd@0
|
50 if clamp_covY, SigmaY = SigmaY; end
|
|
wolffd@0
|
51 if clamp_weights, weightsY = weightsY; end
|
|
wolffd@0
|
52
|
|
wolffd@0
|
53 if nx==1 & ny==1 & regress
|
|
wolffd@0
|
54 P = polyfit(X,Y); % Y = P(1) X^1 + P(2) X^0 = ax + b
|
|
wolffd@0
|
55 assert(approxeq(muY, P(2)))
|
|
wolffd@0
|
56 assert(approxeq(weightsY, P(1)))
|
|
wolffd@0
|
57 end
|
|
wolffd@0
|
58
|
|
wolffd@0
|
59 %%%%%%%% Test
|
|
wolffd@0
|
60 if 0
|
|
wolffd@0
|
61 c1 = randn(2,100); c2 = randn(2,100);
|
|
wolffd@0
|
62 y = c2(1,:); X = [ones(size(c1,2),1) c1'];
|
|
wolffd@0
|
63 b = regress(y(:), X); % stats toolbox
|
|
wolffd@0
|
64 [m,s,w] = linear_regression(c1, y);
|
|
wolffd@0
|
65 assert(approxeq(b(1),m))
|
|
wolffd@0
|
66 assert(approxeq(b(2), w(1)))
|
|
wolffd@0
|
67 assert(approxeq(b(3), w(2)))
|
|
wolffd@0
|
68 end
|
