% Example of least squares fitting
%
% Samuli Siltanen November 2012

% Plot parameters
msize  = 15;
fsize  = 16;
titlefsize = 20;
lwidth = 2;

% Choose amplitude of noise
noiselevel= 3;

% Choose parameters for the line
a = .8;
b = 10;

% Create ideal (y) and noisy (m) data
x = 10:5:100;
y = a*x+b;
m = y+noiselevel*randn(size(y));

% Construct system matrix
A = [x(:),ones(length(x),1)];

% Compute the least squares fit in two ways.

% First: Moore-Penrose pseudoinverse
[U,S,V] = svd(A);
SS      = zeros(size(S.'));
for iii = 1:min(size(S))
    SS(iii,iii) = 1/S(iii,iii);
end
reg1 = V*SS*(U.')*m(:);

% Second: QR decomposition
[Q,R] = qr(A);
m2    = (Q.')*m(:);
reg2  = inv(R(1:size(R,2),:))*m2(1:size(R,2));

% Third: Matlab's simplified use of QR decomposition
% R = triu(qr(A));
% reg2 = R\(R'\(A'*m(:)));
% r = m(:) - A*reg2;
% e = R\(R'\(A'*r));
% reg2 = reg2 + e;


% Plot the result of least squares fit using Moore-Penrose pseudoinverse
figure(2)
clf
t = linspace(10,100,2);
subplot(1,2,1)
plot(t,reg1(1)*t+reg1(2),'b','linewidth',lwidth)
hold on
plot(x,m,'r.','markersize',msize)
axis equal
set(gca,'ytick',[0:20:100],'fontsize',fsize)
set(gca,'xtick',[0:20:100],'fontsize',fsize)
box off
axis([0 110 0 100])
title('Pseudoinverse','fontsize',fsize)



% Plot the result of least squares fit using QR decomposition
subplot(1,2,2)
plot(t,reg1(1)*t+reg1(2),'b','linewidth',lwidth)
hold on
plot(x,m,'r.','markersize',msize)
set(gca,'ytick',[0:20:100],'fontsize',fsize)
set(gca,'yticklabel',{})
set(gca,'xtick',[0:20:100],'fontsize',fsize)
axis equal
axis([0 110 0 100])
box off
title('QR decomposition','fontsize',fsize)


