% Alexandre Trilla (atrilla)
% 22-Jul-2012
% CC BY-NC-SA

% Least Squares regression experiment
% -----------------------------------

% Original line
% y = 0.5x + 1
x = 0:10;
y = 0.5*x + 1;
% Noisy samples
y = y + randn(1, size(x,2))*0.5;
index = ceil(rand*size(x,2));
oulfact = 10;
% Outlier
y(index) = y(index) + oulfact*rand;

disp 'OLS:'

A=[x' ones(size(x,2),1)];
b= y';
AA = A'*A;
Ab = A'*b;
r = pinv(AA)*Ab

plot(x, y, 'k.');

hold on

plot(x, r(1)*x + r(2), 'g');

% ---------------------------

disp 'TLS with mean subtraction:'

mx = mean(x);
my = mean(y);

xm = x - mx;
ym = y - my;
A = [xm' ym']; 
[U,D,V] = svd(A);
r = V(:, size(V,2));
r = [r; (-r(1)*mx - r(2)*my)]

plot(x, (-r(1)/r(2))*x - r(3)/r(2), 'r');

%-------------------------------

disp 'IRWLS:'

rr = robustfit(x',y')

plot(x, rr(2)*x + rr(1), 'b');

legend('Points', 'OLS', 'TLS', 'IRWLS');
title('Least Squares (LS) regression');

hold off
