|
wolffd@0
|
1 function [a, b, error] = weightedRegression(x, z, w)
|
|
wolffd@0
|
2 % [a , b, error] = fitRegression(x, z, w);
|
|
wolffd@0
|
3 % % Weighted scalar linear regression
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % Find a,b to minimize
|
|
wolffd@0
|
6 % error = sum(w * |z - (a*x + b)|^2)
|
|
wolffd@0
|
7 % and x(i) is a scalar
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 if nargin < 3, w = ones(1,length(x)); end
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 w = w(:)';
|
|
wolffd@0
|
12 x = x(:)';
|
|
wolffd@0
|
13 z = z(:)';
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 W = sum(w);
|
|
wolffd@0
|
16 Y = sum(w .* z);
|
|
wolffd@0
|
17 YY = sum(w .* z .* z);
|
|
wolffd@0
|
18 YTY = sum(w .* z .* z);
|
|
wolffd@0
|
19 X = sum(w .* x);
|
|
wolffd@0
|
20 XX = sum(w .* x .* x);
|
|
wolffd@0
|
21 XY = sum(w .* x .* z);
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 [b, a] = clg_Mstep_simple(W, Y, YY, YTY, X, XX, XY);
|
|
wolffd@0
|
24 error = sum(w .* (z - (a*x + b)).^2 );
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 if 0
|
|
wolffd@0
|
27 % demo
|
|
wolffd@0
|
28 seed = 1;
|
|
wolffd@0
|
29 rand('state', seed); randn('state', seed);
|
|
wolffd@0
|
30 x = -10:10;
|
|
wolffd@0
|
31 N = length(x);
|
|
wolffd@0
|
32 noise = randn(1,N);
|
|
wolffd@0
|
33 aTrue = rand(1,1);
|
|
wolffd@0
|
34 bTrue = rand(1,1);
|
|
wolffd@0
|
35 z = aTrue*x + bTrue + noise;
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 w = ones(1,N);
|
|
wolffd@0
|
38 [a, b, err] = weightedRegression(x, z, w);
|
|
wolffd@0
|
39
|
|
wolffd@0
|
40 b2=regress(z(:), [x(:) ones(N,1)]);
|
|
wolffd@0
|
41 assert(approxeq(b,b2(2)))
|
|
wolffd@0
|
42 assert(approxeq(a,b2(1)))
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44 % Make sure we go through x(15) perfectly
|
|
wolffd@0
|
45 w(15) = 1000;
|
|
wolffd@0
|
46 [aW, bW, errW] = weightedRegression(x, z, w);
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48 figure;
|
|
wolffd@0
|
49 plot(x, z, 'ro')
|
|
wolffd@0
|
50 hold on
|
|
wolffd@0
|
51 plot(x, a*x+b, 'bx-')
|
|
wolffd@0
|
52 plot(x, aW*x+bW, 'gs-')
|
|
wolffd@0
|
53 title(sprintf('a=%5.2f, aHat=%5.2f, aWHat=%5.3f, b=%5.2f, bHat=%5.2f, bWHat=%5.3f, err=%5.3f, errW=%5.3f', ...
|
|
wolffd@0
|
54 aTrue, a, aW, bTrue, b, bW, err, errW))
|
|
wolffd@0
|
55 legend('truth', 'ls', 'wls')
|
|
wolffd@0
|
56
|
|
wolffd@0
|
57 end
|
|
wolffd@0
|
58
|
