clear all;
clc;

% dane wejsciowe, P=4 - ilosc obserwacji, N=3 - ilosc zmiennych
X = [2, 1, 3;
    3, 0, 2;
    1, 3, 1;
    4, 4, 0];

P = size(X, 1); %ilosc obserwacji

% transponowanie danych
X = X'; %3x4

% macierz korelacji
R = (X * (X')) / P; %3x3

% wartosci wlasne i wektory wlasne macierzy korelacji - posortowane
% malejaco
[V,D] = eig(R);
[d,I] = sort(diag(D), 'descend');
D = diag(d);
V_sorted = V(:, I); %3x3

K = 2; % ilosc skladnikow glownych

W = (V_sorted (:, 1:K))'; %2x3

% odtworzenie macierzy korelacji
L = diag (d (1:K)); %2x2
R_hat = W' * L * W; %3x2 * 2x2 * 2x3 = 3x3

% transformacja PCA
Y = W * X; % 2x3 * 3x4 = 2x4

% odtworzenie danych wejsciowych
X_hat = W' * Y; % 3x2 * 2x4 = 3x4

% bledy odtworzen
er = norm (R - R_hat) % bezwzgledny
der = er / norm (R) % wzgledny
ex = norm (X - X_hat)
dex = ex / norm (X)