% GaussianDrawsPCA.m
%
% This m-file computes the PCA for n draws from a 2-dim Gaussian PDF.
clear all, close all
% Enter the covariance matrix.
C = [2 1 ; 1 6];
R = chol(C);
N = 1e4;
for i = 1:N
  X(:,i) = R'*randn(2,1);
end
% Compute the PCA of the data corresponding to the columns of X
% Note that svd gives the same decomposition as eig, but the ordering of
% the eigenvalues is in descending, rather than ascending, order, which is
% what we want.
[U,D]=svd(X*X'/N);

% Plot PCA directions
xmin = min(X(1,:)); xmax = max(X(1,:)); 
ymin = min(X(2,:)); ymax = max(X(2,:));
minD = min(xmin,ymin); maxD = max(xmax,ymax);
x = linspace(xmin,xmax); 
phi1 = ( U(2,1)/U(1,1) ) * x; % First PCA vector line
phi2 = ( U(2,2)/U(1,2) ) * x; % Second PCA vector line

figure(1)
 plot(X(1,:),X(2,:),'*')
 hold on, plot(x,phi1,'r'), plot(x,phi2,'g'), hold off
 axis([minD maxD minD maxD])
 legend('data','1st PCA vec','2nd PCA vec')