close all
clear all
%-----------------------------------
 disp('Cork Boring Data: Exe 2.23')
  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'

%Data Matrix Display
cork = cork(:,2:5);
figure(1)
imagesc(cork)

figure(2)
gplotmatrix(cork,[])

figure(3)
andrewsplot(cork)

figure(4)
glyphplot(cork,'glyph','face')


mu = mean(cork)
Sigma = cov(cork)

%Sigma by definition
n=28 %number of observations
J = ones(n,n); 
Sigma1 = 1/(n-1) * (cork' * cork - 1/28 * cork' * J * cork)
H = eye(n) - 1/n * J; %
Sigma2 =   1/(n-1) * cork' * H * cork 
% Sigma, Sigma1 and Sigma2 are the same.

D = diag( sqrt(diag(Sigma))) ;  %diag(Sigma) is a vector, 
%repeated diag  is needed for diagonal matrix.
%MatrixD connects correlation (R)  and covariance (Sigma)
%    matrices as Sigma = D R D.
% STANDARDIZED DATA SET
%
% Center cork first...
corkc = cork -  repmat( mean(cork),n,1 );
Y = (inv(D) * corkc' )'
% or equivalently Y = corkc * inv(D)
cov(Y)  %is equal to R=corr(Y).

% DATA SET WITH DECORRELATED VARIABLES
%Mahalanobis Transform 
M = sqrtm( inv(Sigma));

Z = corkc * M;

cov(Z)

figure(5)
subplot(1,3,1)
imagesc(cork)
subplot(1,3,2)
imagesc(Y)
subplot(1,3,3)
imagesc(Z)

%%
%PRINCIPAL COMPONENTS
close all force
%Rao's Cork Data components (p=4) are similar in scale.
% Therefore Cov matrix is used.
[n p]=size(cork);

% covariance  matrix
C=cov(cork)  %  
% MATLAB built in. It outputs  PC's (coeffs), transformed data (scores), and
% variances (latent). 
[coeffs, scores, latent] = pca(cork);

 %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 cork*V
[V D]=eigs(C,4);
coeffs2=V ; 
latent2=diag(D);   
scores2= cork*V; 



%plot
figure;
pareto(latent)

latent(1)/sum(latent)  % 90.62% variance in 1st PC

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')




