wolffd@0: function [a, b, error] = weightedRegression(x, z, w)
wolffd@0: % [a , b, error] = fitRegression(x, z, w);
wolffd@0: % % Weighted scalar linear regression
wolffd@0: %
wolffd@0: % Find a,b to minimize
wolffd@0: % error = sum(w * |z - (a*x + b)|^2) 
wolffd@0: % and x(i) is a scalar
wolffd@0: 
wolffd@0: if nargin < 3, w = ones(1,length(x)); end
wolffd@0: 
wolffd@0: w = w(:)';
wolffd@0: x = x(:)';
wolffd@0: z = z(:)';
wolffd@0: 
wolffd@0: W = sum(w);
wolffd@0: Y = sum(w .* z);
wolffd@0: YY = sum(w .* z .* z);
wolffd@0: YTY = sum(w .* z .* z);
wolffd@0: X = sum(w .* x);
wolffd@0: XX = sum(w .* x .* x);
wolffd@0: XY = sum(w .* x .* z);
wolffd@0: 
wolffd@0: [b, a] = clg_Mstep_simple(W, Y, YY, YTY, X, XX, XY);
wolffd@0: error = sum(w .* (z - (a*x + b)).^2 );
wolffd@0: 
wolffd@0: if 0
wolffd@0:   % demo
wolffd@0:   seed = 1;
wolffd@0:   rand('state', seed);   randn('state', seed);
wolffd@0:   x = -10:10;
wolffd@0:   N = length(x);
wolffd@0:   noise = randn(1,N);
wolffd@0:   aTrue = rand(1,1);
wolffd@0:   bTrue = rand(1,1);
wolffd@0:   z = aTrue*x + bTrue + noise;
wolffd@0:   
wolffd@0:   w = ones(1,N);
wolffd@0:   [a, b, err] = weightedRegression(x, z, w);
wolffd@0:   
wolffd@0:   b2=regress(z(:), [x(:) ones(N,1)]);
wolffd@0:   assert(approxeq(b,b2(2)))
wolffd@0:   assert(approxeq(a,b2(1)))
wolffd@0: 
wolffd@0:   % Make sure we go through x(15) perfectly
wolffd@0:   w(15) = 1000;
wolffd@0:   [aW, bW, errW] = weightedRegression(x, z, w);
wolffd@0: 
wolffd@0:   figure;
wolffd@0:   plot(x, z, 'ro')
wolffd@0:   hold on
wolffd@0:   plot(x, a*x+b, 'bx-')
wolffd@0:   plot(x, aW*x+bW, 'gs-')
wolffd@0:   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: 		aTrue, a, aW, bTrue, b, bW, err, errW))
wolffd@0:   legend('truth', 'ls', 'wls')
wolffd@0:   
wolffd@0: end
wolffd@0: