% pca_ukazka - skript na ukazku analyzy hlavnich komponent (PCA)
%
% Autor: Eva Koritakova (Janousova)
% Masarykova univerzita, Brno
% 14/04/2015, modified 07/11/2017
%
% ----------------------------------------------------------------------- %

%nacteme datovy soubor (1. list z dokumentu Microsoft Excel)
[num, txt, raw] = xlsread('Data_neuro.xlsx',1);
data = num(:,24:29); % vyber 6 promennych s objemy mozkovych struktur
clear num txt raw


%% PCA:
[coeff,score,latent] = pca(data);
% coeff - matice vlastnich vektoru (vl. vektory jsou ve sloupcich)
% score - souradnice subjektu v novem prostoru (tzn. v radcich jsou
%        subjekty, ve sloupcich hlavni komponenty serazene podle vl. cisel)
% latent - vlastni cisla

%Pokud bychom si chteli souradnice pomoci vl. vektoru spocitat sami:
score2=data*coeff; % vynasobeni dat matici vlastnich vektoru
score2=score2-ones(size(data,1),1)*mean(score2); % centrovani (tzn. odecteni prumeru z kazde hlavni komponenty)

                                      
%% Vypocet PCA pomoci funkce eig:

% 1. Vypocet kovariancni matice:
kov_mat=cov(data);

% 2. Vypocet vlastnich cisel a vlastnich vektoru matice:
[V,lambda]=eig(kov_mat); % V...matice ortonormalnich vlastnich vektoru,lambda...matice,ktera ma na diagonale vl. cisla
lambda=diag(lambda); % naskladani vlastnich cisel do sloupcoveho vektoru

% 3. Vyber komponent a vytvoreni feature vektoru (vektoru vlastnich vektoru):
lambda=lambda(end:-1:1,1); % serazeni vl.cisel od nejvetsiho po nejmensi (z fce eig vypadnou serazene od nejmensiho po nejvetsi)
V=V(:,end:-1:1); % serazeni vlastnich vektoru podle serazeni vlastnich cisel (matice V je totozna jako matice coeff vypocitana pomoci funkce pca)

% 4. Vypocet procent vycerpane variability a vykresleni sutinoveho grafu:
sum_lambda=cumsum(lambda)/sum(lambda); % normalizace vektoru kumulovanych souctu vlastnich cisel
figure, plot(lambda) % sutinovy graf
figure, plot(sum_lambda) % zobrazeni kumulativniho souctu vlastnich cisel,aby urceni poctu komponent bylo prehlednejsi (vidime snaze, kolik procent variability bude zachovano)
ylim([0 1])

% 7. Vytvoreni nove datove sady (eigenimages):
dataY=data*V; % pouzijeme pro transformaci dat vsechny hl. komponenty
dataY=dataY-ones(size(data,1),1)*mean(dataY); % centrovani (tzn. odecteni prumeru z kazde hlavni komponenty)
% matice dataY je totozna s matici score a score2

% 8. Korelace s puvodnimi promennymi:
korelace=corr(data,dataY);

% kdybychom chteli korelace pocitat postupne, provadeli bychom to takto:
corr(data(:,1),dataY(:,1)) % korelace objemu hipokampu s 1. hl. komponentou - verze 1
    (V(1,1)*sqrt(lambda(1)))/sqrt(kov_mat(1,1)) % korelace objemu hipokampu s 1. hl. komponentou - verze 2
corr(data(:,2),dataY(:,1)) % korelace objemu amygdaly s 1. hl. komponentou - verze 1
    (V(2,1)*sqrt(lambda(1)))/sqrt(kov_mat(2,2)) % verze 2 - je to V(2,1) a ne V(1,2) proto, ze vl. vektory jsou ve sloupcich
    


