% Work through a simple example of principal component analysis (PCA)
% without explicit use of svd
%
% Samuli Siltanen November 2012

% Plot parameters
msize  = 15;
lwidth = 1.5;
fsize  = 16;

% Construct the data
xvec = [2.4 0.6 2.1 1.7 3.2 2.4 2.1 .8 1.6 1.2];
yvec = [2.5 0.6 2.8 2.3 3.1 2.6 1.5 1.2 1.5 1];

% Subtract means
xvec0 = xvec-mean(xvec);
yvec0 = yvec-mean(yvec);

% Construct the data matrix where the data vectors appear as rows
XT = [xvec0(:),yvec0(:)];

% Calculate covariance matrix of the data
covmat = cov(XT);

% Compute eigenvalues and eigenvectors of the covariance matrix.
% We change the order of eigenvalues so that the bigger comes first.
[P,D] = eig(covmat);
P = [P(:,2),P(:,1)];
D = diag(flipud(diag(D)));

% Transform the data
dmat2 = inv(P)*(XT.');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Show the data

% Determine axis limits 
MINx = min([min(xvec),min(xvec0),min(dmat2(1,:))]);
MAXx = max([max(xvec),max(xvec0),max(dmat2(1,:))]);
MINy = min([min(yvec),min(yvec0),min(dmat2(2,:))]);
MAXy = max([max(yvec),max(yvec0),max(dmat2(2,:))]);

figure(1)
clf
subplot(1,3,1)
plot(xvec,yvec,'b.','markersize',msize)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])
title(['Original data'],'fontsize',fsize)

pause

subplot(1,3,2)
plot(xvec0,yvec0,'b.','markersize',msize)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])
set(gca,'yticklabel',{})
title(['Zero-mean data'],'fontsize',fsize)

pause

% Show principal directions
subplot(1,3,2)
plot(xvec0,yvec0,'b.','markersize',msize)
hold on
plot([0 2*sqrt(D(1,1))*P(1,1)],[0 2*sqrt(D(1,1))*P(2,1)],'r','linewidth',lwidth)
plot([0 2*sqrt(D(2,2))*P(1,2)],[0 2*sqrt(D(2,2))*P(2,2)],'r','linewidth',lwidth)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])
set(gca,'yticklabel',{})
title(['Zero-mean data'],'fontsize',fsize)

pause

subplot(1,3,3)
plot(dmat2(1,:),dmat2(2,:),'b.','markersize',msize)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])
set(gca,'yticklabel',{})
title(['Transformed data'],'fontsize',fsize)

