PCA Lab

clear;
close all;
load("iris.mat");
setosa_label = 1;
virginica_label = 2;
versicolor_label = 3;

setosa = [ones(size(setosa,1),1)*setosa_label setosa];
virginica = [ones(size(virginica,1),1)*virginica_label virginica];
versicolor = [ones(size(versicolor,1),1)*versicolor_label versicolor];
data = [setosa ; virginica ; versicolor];

function [X] = centerAndScale(data)
  means = mean(data);
  stdv = std(data);
  X = zeros(size(data));
  for i=1:size(data,2)
    for j=1:size(data,1)
        X(j,i) = ( means(:,i) - data(j,i))/stdv(:,i);
    end
  end
end

data(:,2:end) = centerAndScale(data(:,2:end));

mean_orig_1 = mean(data(1:50,2:end))
mean_orig_2 = mean(data(51:100,2:end))
mean_orig_3 = mean(data(101:150,2:end))

[U, S, V] = svd(data(:,2:end),0)
Ur = U * S;
data(:,2:end)=Ur;

mean_1 = mean(data(1:50,2:end))
mean_2 = mean(data(51:100,2:end))
mean_3 = mean(data(101:150,2:end))

Ur_1_cov = cov(data(1:50,2:end))
Ur_2_cov = cov(data(51:100,2:end))
Ur_3_cov = cov(data(101:150,2:end))

function [weights, cumulative_weights] = scree(S)
  S2 = S^2;
  weights = zeros(size(S,2),1);
  sumS2 = sum(sum(S2));
  weightsum = 0;

  for i=1:size(S,2)
    weights(i) = S2(i,i)/sumS2;
    weightsum = weightsum + weights(i);
    cumulative_weights(i) = weightsum;
  end
end

[weights, cumulative_weights] = scree(S)
figure;
plot(weights,'x:b', "markersize", 10, "linewidth", 5, "linestyle", '-');
set(gca, "fontsize", 28)
xlabel("Principal Component", "Fontsize",36, "Fontweight","Bold");
ylabel("Proportion Of Variance", "Fontsize",36, "Fontweight","Bold");
title('Scree Plot',"Fontsize",48, "Fontweight","Bold");
grid;
figure;
plot(cumulative_weights,'x:b', "markersize", 10, "linewidth", 5, "linestyle", '-');
set(gca, "fontsize", 28)
grid;
xlabel("Principal Component", "Fontsize",36, "Fontweight","Bold");
ylabel("Cumulative Proportion Of Variance", "Fontsize",36, "Fontweight","Bold");
title('Cumulative Scree Plot',"Fontsize",48, "Fontweight","Bold");

setosa = data(1:50,:);
setosa_mean = mean(setosa(:,2:end));
setosa_covariance = cov(setosa(:,2:end));
virginica = data(51:100,:);
virginica_mean = mean(virginica(:,2:end));
virginica_covariance = cov(virginica(:,2:end));
versicolor = data(101:150,:);
versicolor_mean = (mean(versicolor(:,2:end)));
versicolor_covariance = (cov(versicolor(:,2:end)));

function [Vsquare] = squareSign(V)
  for i=1:size(V,2)
    for j=1:size(V,2)
        Vsquare(i,j) = V(i,j)^2;
        if V(i,j)<0
            Vsquare(i,j) = Vsquare(i,j)*-1;
        else
            Vsquare(i,j) = Vsquare(i,j)*1;
        end
    end
  end
end

function [] = plotLoading(Vsquare, idx)
  figure;
  bar(Vsquare(:,idx),0.5);
  grid;
  set(gca, "fontsize", 28)
  ymin = min(Vsquare(:,idx)) + (min(Vsquare(:,idx))/10);
  ymax = max(Vsquare(:,idx)) + (max(Vsquare(:,idx))/10);
  axis([0 size(Vsquare, 2)+1 ymin ymax]);

  xlabel('Feature index', "fontsize", 36);

  ylabel('Importance of feature', "fontsize", 36);

  [chart_title, ERRMSG] = sprintf('Loading Vector %d',idx);

  title(chart_title, "fontsize", 48);
end

[Vsquare] = squareSign(V)
plotLoading(Vsquare, 1)
plotLoading(Vsquare, 2)
plotLoading(Vsquare, 3)
plotLoading(Vsquare, 4)

function [] = plotScatter3(setosa, versicolor, virginica, idx, idy, idz, labels)
  figure;

  scatter3(setosa(:,idx), setosa(:,idy), setosa(:,idz), 5, 'r')
  hold on;
  scatter3(versicolor(:,idx), versicolor(:,idy), versicolor(:,idz), 5, 'g')
  scatter3(virginica(:,idx), virginica(:,idy), virginica(:,idz), 5, 'b')
  set(gca, "fontsize", 28)
  xlabel(labels{idx - 1}, "fontsize", 36)
  ylabel(labels{idy - 1}, "fontsize", 36)
  zlabel(labels{idz - 1}, "fontsize", 36)
  title([labels{idx - 1} " vs. " labels{idy - 1} " vs. " labels{idz - 1}], "fontsize", 48)

  hold off;
end

# Plot original data
original_labels = cellstr(['Sepal Width'; 'Sepal Length'; 'Petal Width'; 'Petal Length']);

plotScatter3(setosa, versicolor, virginica, 2, 3, 4, original_labels);
plotScatter3(setosa, versicolor, virginica, 2, 3, 5, original_labels);
plotScatter3(setosa, versicolor, virginica, 2, 4, 5, original_labels);
plotScatter3(setosa, versicolor, virginica, 3, 4, 5, original_labels);

data_setosa = data(1:50, :);
data_versicolor = data(51:100, :);
data_virginica = data(101:150, :);

# Plot PCA data
pca_labels = cellstr(['PC1'; 'PC2'; 'PC3'; 'PC4']);

plotScatter3(data_setosa, data_versicolor, data_virginica, 2, 3, 4, pca_labels);
plotScatter3(data_setosa, data_versicolor, data_virginica, 2, 3, 5, pca_labels);
plotScatter3(data_setosa, data_versicolor, data_virginica, 2, 4, 5, pca_labels);
plotScatter3(data_setosa, data_versicolor, data_virginica, 3, 4, 5, pca_labels);