%
%
% Linear least squares demo
%

% Load the Filip data set
% (see http://www.itl.nist.gov/div898/strd/lls/data/Filip.shtml)
load Filip.dat
y = Filip(:,1);
x = Filip(:,2);


% perturb measurements
%y = y.*(1+0.01*randn(size(y)));


% Plot the data
plot(x,y,'.'); hold on
title(' NIST Filip Data Set')
xlabel('x');
ylabel('y')

% Fit polynomial of degree 10, i.e. find beta_0, ..., beta_10
% such that  
%   y approx beta_0 + beta_1 x + ... + beta_10 x^10
% in the least squares sense

degree = 10;

% Set up least squares problem
b = y;
A = ones(size(y,1),degree+1);
for i = 2:degree+1
    A(:,i) = A(:,i-1).*x;
end


% Solve using Normal equation approach:
AtA   = A'*A;
[R,p] = chol(AtA);
if( p > 0 )
    disp(['p = ',int2str(p),' > 0 '])
    disp('Cholesky decomposition could not be performed, use backslash')
    beta = AtA \ (A'*b);
else
    beta = R\(R'\(A'*b));
end

disp('polynomial coefficients computed using the normal equation approach')
disp(beta)
%plot the polynomial
figure(1)
xx = [-9:0.1:-3]';
x2 = ones(size(xx));
yy = beta(1)*x2;
for i = 1:degree
    x2 = x2.*xx;
    yy = yy + beta(i+1)*x2;
end
plot(xx,yy,'-r')

disp('Hit RETURN to continue....'); pause



% Solve using QR decomposition:
[Q,R,P] = qr(A);
figure(2)
semilogy(abs(diag(R)),'*')
title('absolute values of diagonals of R')
d       = Q'*b;
beta    = P*(R(1:degree+1,1:degree+1) \ d(1:degree+1));
disp('polynomial coefficients computed using the QR decomposition')
disp(beta)

%plot the polynomial
figure(1)
xx = [-9:0.1:-3]';
x2 = ones(size(xx));
yy = beta(1)*x2;
for i = 1:degree
    x2 = x2.*xx;
    yy = yy + beta(i+1)*x2;
end
plot(xx,yy,'-g');


disp('Hit RETURN to continue....'); pause

% Solve using SVD:
[U,S,V] = svd(A);
s       = diag(S);
figure(2)
semilogy(s,'*')
title('singular values of A')
figure(3)
plot(V)
title('singular vectos of A')
d       = U'*b;
beta    = V* (d(1:degree+1)./s);
disp('polynomial coefficients computed using the SVD')
disp(beta)

%plot the polynomial
figure(1)
xx = [-9:0.1:-3]';
x2 = ones(size(xx));
yy = beta(1)*x2;
for i = 1:degree
    x2 = x2.*xx;
    yy = yy + beta(i+1)*x2;
end
plot(xx,yy,'--m'); hold off
legend('data','polynomial computed using NE',...
       'polynomial computed using QR','polynomial computed using SVD')