% 
% problem epage: 69
% example 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

x1 = randn(50,1)*100; 
x2 = randn(50,2); 
X  = [x1,x2]; 

n = 50; p = 3; 

%--
% perform PCA on the covariance matrix: 
%--
covm = cov(X); 

[eigvec,eigval] = eig(covm); 
eigval = diag(eigval); 

eigval = flipud(eigval);   % Order in decending order: 
eigvec = eigvec(:,p:-1:1); % permute the order of the columns to match: 

% display the eigenvalues:
fprintf('%20.6f',eigval); fprintf('\n'); 

% lets project onto the dimension with the largest variance:
v1 = X * eigvec(:,1);

% lets see how close this is to x1 (modulo) sign: 
[ -X(1:5,1), v1(1:5) ]  % print a comparison of the two data points 
fprintf('distance to x1 = %10.6f\n', norm( -X(:,1)-v1 ) ); 

% Do a scree plot: 
fh=figure; 
ph1=plot(1:length(eigval),eigval,'ko-'); hold on; 
title( 'cov scree plot' ); grid on; 

%--
% perform PCA on the correlation matrix: 
%--
corm = corrcoef(X); 

[eigvec,eigval] = eig(corm); 
eigval = diag(eigval); 

eigval = flipud(eigval);   % Order in decending order: 
eigvec = eigvec(:,p:-1:1); % permute the order of the columns to match: 

% display the values 
fprintf('%20.6f',eigval); fprintf('\n'); 

% lets project onto the dimension with the largest variance:
v1 = X * eigvec(:,1);

% lets see how close this is to x1 (modulo) sign: 
fprintf('distance to x1 = %10.6f\n', norm( X(:,1)-v1 )); 

% Do a scree plot: 
figure; 
ph2=plot(1:length(eigval),eigval,'ko-'); 
%legend( [ph1,ph2], {'PCA cov', 'PCA corr'}); 
title( 'correlation based scree plot' );