% 
% Epage 70 
% 
% 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.
% 
%-----

addpath( '../../Code/eda_data' ); 

close all; drawnow; 
%clc; 
clear; 

% Pt a: 
n=30; p=3;
%n=30000; p=3;
dataOrig = randn(n,p); 
datac = dataOrig - repmat(sum(dataOrig)/n,n,1); 
% Find the covariance matrix:
covm = cov(datac); 
[eigvec,eigval] = eig(covm); 
eigval = diag(eigval); 
% Order in decending order: 
eigval = flipud(eigval); 
eigvec = eigvec(:,p:-1:1); 

fprintf('total of all eigenvalues = %10.6f\n', sum( eigval ) ); 
fprintf('trace of covm = %10.6f\n', sum( diag( covm ) ) ); 

ZpcScore_pta = datac * eigvec; 
fprintf( 'mean of the projected PCA scores Part (a) = \n' ); 
mean( ZpcScore_pta ) 

% Pt b: 
data = dataOrig + repmat( [ 2 2 2 ], [size(dataOrig,1),1] ); 
covm = cov(data); 
[eigvec,eigval] = eig(covm); 
eigval = diag(eigval); 
% Order in decending order: 
eigval = flipud(eigval); 
eigvec = eigvec(:,p:-1:1); 

ZpcScore_ptb = data * eigvec; 
fprintf( 'mean of the projected PCA scores Part (b) = \n' ); 
mean( ZpcScore_ptb ) 
fprintf( 'the mean [2,2,2] projected on the PCA scores = \n' ); 
[ 2, 2, 2 ] * eigvec 

% Pt c: 
% datac is already centered, scale by the stds: 
s = std( datac ); 
datacp = datac ./ repmat( s, size(datac,1), 1 ); 
covt = cov( datacp ); 

% Compare the covariance/correlation: 
fprintf( 'the difference between covariance and correlation = %10.6f\n', ...
         norm( covt - corrcoef( dataOrig ) ) ); 

% PT d:
ZpcScore = data * eigvec; 
fprintf( 'the covariance of the PC scores is given by \n' );
cov( ZpcScore )