% Work through a simple example of principal component analysis (PCA)
% using singular value decomposition (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
m1    = mean(xvec);
m2    = mean(yvec);
xvec0 = xvec-m1;
yvec0 = yvec-m2;

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

% Calculate SVD
X       = XT.';   % Taking the transpose
[W,S,V] = svd(X); % This results in X=W*S*V.'

% Transform the data using the formula YT = XT*W = V*ST
YT = V*(S.');
Y  = YT.';

% Dimension reduction
SS      = S;
SS(2,2) = 0;
Y2T     = V*(SS.');
Y2      = Y2T.';

% Transform back
X2 = W*Y2;

% Add back the mean value
n  = length(xvec);
X3 = X2+[m1*ones(1,n);m2*ones(1,n)];

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

% Show the data

% Determine axis limits 
MINx = min([min(xvec),min(xvec0),min(Y(1,:))]);
MAXx = max([max(xvec),max(xvec0),max(Y(1,:))]);
MINy = min([min(yvec),min(yvec0),min(Y(2,:))]);
MAXy = max([max(yvec),max(yvec0),max(Y(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(S(1,1))*W(1,1)],[0 2*sqrt(D(1,1))*W(2,1)],'r','linewidth',lwidth)
%plot([0 2*sqrt(S(2,2))*W(1,2)],[0 2*sqrt(D(2,2))*W(2,2)],'r','linewidth',lwidth)
plot([0 S(1,1)*W(1,1)],[0 S(1,1)*W(2,1)],'r','linewidth',lwidth)
plot([0 S(2,2)*W(1,2)],[0 S(2,2)*W(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(Y(1,:),Y(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)

pause

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

pause

subplot(1,3,2)
plot(xvec0,yvec0,'b.','markersize',msize)
hold on
plot(X2(1,:),X2(2,:),'k.','markersize',msize)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])
set(gca,'yticklabel',{})

pause

subplot(1,3,1)
plot(xvec,yvec,'b.','markersize',msize)
hold on
plot(X3(1,:),X3(2,:),'k.','markersize',msize)
axis equal
axis([1.1*MINx,1.1*MAXx,1.1*MINy,1.1*MAXy])

