X = load('Grade.csv');    
X = X';                    % rows of X are now grades on particular Final problems, for n = 113 students 
[p,n] = size(X);   % number of vars

%============================ Part I :  using Covariance matrix approach

figure(1)
Ca = cov(X');         % one way to find the covariance matrix            
   subplot(1,2,1)
   imagesc (Ca);  colormap(gray)
   title('covariance matrix')

mx = mean(X,2);
X = X - mx*ones(1,n);           % mean-corrected X 
C = X*X';             % another way to find covariance matrix, up to a constant 
  subplot(1,2,2)
  imagesc(C);      colormap(gray)

Corr = corrcoef(X');
Corr(1:3,1:3)
 imagesc(Corr);
    title('correlation matrix')

[U,D] = eig(C);

figure(2)
plot((p:-1:1)',diag(D)); 
   xlabel('PC number')
   ylabel('Eigenvalue size')
  
PC1 = U(:,p);        % Matlab lists the eigenvalues in *increasing* order
PC2 = U(:,p-1);
PC3 = U(:,p-2);

figure(3)
plot(PC1, PC2, 'w*');  bookfonts;
  xlabel('PC1'); ylabel('PC2');  title('Contributions from questions 1-21')
   hold on
   for i = 1:p,
     text(PC1(i), PC2(i),int2str(i),...
         'Fontsize',18)  % shows the problem number 
   end
   hold off

figure(4) 
plot(PC2, PC3,'w*'); bookfonts;
   hold on
   for i = 1:p,
     text(PC2(i), PC3(i),int2str(i),...
         'Fontsize',18)  % shows the problem number 
   end
   hold off


%% now compute the PC1, PC2 for all the students in X 

 A  = U'*X;      % this way we define a "score" matrix
                  % whose rows are the scores for each student on d PC's 
figure(5)                  
plot(A(p,:), A(p-1,:),'*')   % this the plot of first and second PC's  
                    % not an ellipse because many students got close to
                    % 100% -- see the narrowing on the right 
   title('Students'' PC1 and PC2 scores'); bookfonts                   
   xlabel('PC1 score'); ylabel('PC2 score')                    
   % students who did badly mostly suffered on the advanced questions: 
   %   their PC2 values are negative
   
  
  
%% Part II  use of SVD decomposition instead of correlation/eigenvalue approach

[U0,S,V] = svd(X,0);  

 figure(6)
 plot(diag(S))     % compare this with eigenvalues D, Fig. 2

 PC1 = U0(:,1);
 PC2 = U0(:,2);
 
 plot(PC1, PC2,'*')     % compare with Fig. 3
 
 A = S*V'; 
 A1 = A(1,:);  A2 = A(2,:);   % compare with Fig. 5
 plot(A1,A2,'*')
 
 A3 = A(3,:); 
 plot3(A1,A2,A3,'*')