close all
clear all
%-----------------------------------
 disp('Cork Boring Data PCA')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 17);
  fs = 17;
  msize = 5;
%  Cork Boring Data    
%  The cork oak (Quercus suber L), a noble
%  tree with very special characteristics that grows in Mediterranean.
%  This natural tissue, has unique qualities: light weight, elasticity,
%  insulation and impermeability, fire retardancy, resistance to abrasion, etc.
%  The data measuring cork boring of tress given in Rao (1948) consist  of the
%  weights (in centigrams) of cork boring in four directions
% (north, east, south, and west) for 28 trees.
 
load 'C:\STAT\Descriptive\Descriptivedat\cork.dat'

close all force
%Rao's Cork Data components (p=4) are similar in scale-> use cov matrix
[n p]=size(cork);
% covariance  matrix
C=cov(cork)  %  
% MATLAB built in. Gives PC (coeffs), transformed data (scores), and
% variances (latent). 

[coeffs, scores, latent] = pca(cork);
coe1=coeffs(:,1)  %first PC
var1= latent         %Variances split among the PC's
sco1= scores(:,1) %1st coordinate of the data in PC domain. Best univariate summary of the data
                             %in the sense of compressing variance.

 %Equivalent task using eig, eigs. Note that eigenvectors (columns of V) are PCs,
 % and variances are on the diagonal of D.
 % Scores are recovered as Z*V
[V D]=eigs(C,4);
coe2=V(:,1) %V=coeffs
var2=diag(D)  %latent
C22= cork*V;  %scores
sco2=C22(:,1)

%Equivalent task using the most rudimentary operations via the
%singular value decomposition of data matrix.   
%cork = U L  V,  D=1/(n-1)*L^2, V is the matrix of eigenvectors.

cork=cork - repmat(mean(cork), n,1);
[u lam v]=svd(cork);
coe3= v(:,1)  %v=coeffs
var2= diag(lam.^2/(n-1))  %latent
C33= cork * v;  %scores
sco3=C33(:,1)


%plot
figure;
pareto(var1)
var1(1)/sum(var1)
figure;
subplot(211) 
plot(cork(:,1),cork(:,2),'ko')
xlabel('1st coord')
ylabel('2nd coord')
subplot(212)
plot(scores(:,1), scores(:,2),'ko')
xlabel('1st PC')
ylabel('2nd PC')




