% pca_demo.m
% Explain the Principal Component Aanlysis
% Dr. David Pan, UAH

close all
clear all

N = 100;


m1 = [3, 3]';   % Mean vector
cov1 = [2 1; 1 2]; % Covariance matrix
rng default  % For reproducibility
Xraw = mvnrnd(m1,cov1,N);

% Export the dataset to a text file (to be imported by sklearn)
writematrix(Xraw, 'pca.csv');


% Center the data
mu = mean(Xraw);
X = Xraw - mu;
mean(X)

figure;
%plot(Xraw(:,1),Xraw(:,2),'b.','MarkerSize',12)
plot(X(:,1),X(:,2),'b.','MarkerSize',12)
grid

% Empirical sample covariance matrix
(X'*X)/(N-1)

% Same as using cov()
C = cov(X)

% Find out the eigenvalues and eigenvectors
[V,D] = eig(C);
C*V
V*D

V1 = V(:,1);
V2 = V(:,2);
% Display the eigenvectors for data projection
hor = min(X(:,1)): 0.01: max(X(:,1));
ver = V1(2)/V1(1)*hor;
plot(hor, ver, 'g');
axis equal

hold on;
hor = min(X(:,1)): 0.01: max(X(:,1));
ver = V2(2)/V2(1)*hor;
plot(hor, ver, 'r');


% Variance after projection
Y1 = X*V1;
var(Y1)
figure;  scatter(Y1,zeros(100,1))

Y2 = X*V2;
var(Y2)
figure;  scatter(Y2,zeros(100,1))
grid

% Correlation between projections


% Using SVD
[A,S,B]= svd(X/sqrt(N-1),'econ');
S
% Singular values are the same as the square root of the eigenvalues of C
sqrt(D)
% Compare the right singular vector B with the eigenvectors of C

% Use pca function
[coeff,score,latent] = pca(Xraw);

% Reconstruction using only the principal components
Y1_drop = zeros(N,1);  % Setting the first component to zeros
Y = [Y1_drop, Y2];
%Y = [Y1, Y2];  % Keep both components
% inv(V) = V = V'
X_rec = Y*V;

figure;
plot(X_rec(:,1),X_rec(:,2),'b.','MarkerSize',12)
grid

diff = (X_rec - X);
diff_sq = diff(:,1).^2 + diff(:,2).^2;
% Average distortion (mean square error)
sum(diff_sq)/(N-1)   % Same as the eigenvalue of the 1st component

% Repeat the above steps using the score matrix
score_truncated = score;
score_truncated(:,2) = 0;
X_rec2 = score_truncated * coeff';
figure;
plot(X(:,1),X(:,2),'b.','MarkerSize',12)
hold on;
plot(X_rec2(:,1),X_rec2(:,2),'r.','MarkerSize',12)
legend('Original Data', 'Reconstructed Data','Location','SE');
grid
axis equal

diff = (X_rec2 - X);
diff_sq = diff(:,1).^2 + diff(:,2).^2;
% Average distortion (mean square error)
sum(diff_sq)/(N-1) 