clear all
close all


 disp('Multivariate Data: Body Fat Example')
 lw = 3; 
 set(0, 'DefaultAxesFontSize', 16);
 fs = 15;
 msize = 10;


load('C:\STAT\Descriptive\Descriptivedat\fat.dat')
 % load('fat.dat')
casen = fat(:,1);
broz = fat(:,2);
siri = fat(:,3); 
densi = fat(:,4); 
age = fat(:,5);
weight = fat(:,6);   
height = fat(:,7);   
adiposi = fat(:,8); 
ffwei = fat(:,9);   
neck = fat(:,10); 
chest = fat(:,11);
abdomen  = fat(:,12);   
hip = fat(:,13); 
thigh = fat(:,14);
knee = fat(:,15);   
ankle = fat(:,16);
biceps = fat(:,17);
forearm = fat(:,18);     
wrist = fat(:,19);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = [broz  densi weight adiposi biceps];

varNames = {'broz'; 'densi'; 'weight'; 'adiposi'; 'biceps'};

% covariances and correlations

%%
% mean, covariance and correlation
%              matrices
Xbar = mean(X)
S = cov(X)
R = corr(X)

%%
% S, R  by ``hand''
[n p]=size(X);
H = eye(n) - 1/n * ones(n,1)*ones(1,n);
S = 1/(n-1) * X' * H * X;
stds = sqrt(diag(S));
D = diag(stds);
R = inv(D) * S * inv(D) 
%S and R here coincide with S and R 
% calculated by built in finctions cov and cor.

%%

Xc= X- repmat(mean(X),n,1);  %center X
% subtract component means from variables in each observation.

% there are two interesting transformations
Y =  Xc * inv(D); % the scaling transform for 
% which S=R
cov(Y) %is corr(X)

% and Mahalanobis transform 
M = sqrtm(inv(S))  %sqrtm is a sqare root of matrix

Z = Xc * M; %Z has uncorrelated components

cov(Z) %should be the identity matrix
%%
figure(1)
imagesc(X)
colorbar
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\fat5X.eps'

figure(2)
imagesc(Y)
colorbar
%print -depsc 'C:\Springer\Descriptive\Descriptiveeps\fat5Y.eps'


figure(3)
imagesc(Z)
colorbar
% print -depsc 'C:\Springer\Descriptive\Descriptiveeps\fat5Z.eps'

%%
%Principal Components 
close all force
%Take all non-invasive body measures
clear X %clear old X
X=[age weight height adiposi neck chest abdomen hip thigh knee ankle biceps forearm wrist]
[n p]=size(X);

Y=zscore(X); % components of matrix X are on very different scales, so ZSCORE matrix is used.
% correlation  matrix
C=cov(Y)  % cov(Y) is the same as corr(X)

% MATLAB built in. Gives PC (coeffs), transformed data (scores), and
% variances (latent). When Y is used -- sum(variances)=p=14.
[coeffs, scores, latent] = pca(Y);
A1=coeffs(:,1)  %first PC
B1= latent         %Vars
C1= 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 Y*V
[V D]=eigs(C,14);
A2=V(:,1)
B2=diag(D)
C22= Y*V; 
C2=C22(:,1)


%Equivalent task using the most primitive singular value decomposition
%of data.  Y= U S V,  D=S^2, V is matrix of eigenvectors.

[u s v]=svd(Y,'econ');
A3= v(:,1)
B3= diag(s.^2/(n-1))
C33= Y * v;  
C3=C33(:,1)

% A1, A2,and  A3 coincide;
% B1, B2, and B3 coincide;  and
% C1, C2, and C3 coincide

figure;
scatter(C2, broz)
xlabel('FirstPrincipal Component')
ylabel('Brozek Index')

 


