# Studijny program: "Pravdepodobnost a matematicka statistika"
# Vyucujuci: Radoslav Harman, KAMS FMFI UK Bratislava
#
# Metoda hlavnych komponentov (Principal components analysis, PCA)

# Na predmete AZKD pouzivame PCA najma na vizualizaciu a explorativnu analyzu
# PCA je mozne pouzit aj na redukciu poctu premennych

### Najprv si zopakujme spektralny rozklad symetrickej matice pomocou R-ka
A <- matrix(rnorm(25), ncol = 5)
A <- A + t(A)
res <- eigen(A)

lam <- res$values; lam
U <- res$vectors; U

# Sedi to s teoriou? Vsimnime si tiez numericke odchylky
U %*% t(U)
A - U %*% diag(lam) %*% t(U)

# Analyzujme data decathlon 2012 z kniznice scar
library(scar)
data(decathlon)
dec <- decathlon; dec
attach(dec) 

# Pozrime sa na (ne)normalitu jednej premennej
qqnorm(LJ)
shapiro.test(LJ)

# Pozrime sa na (ne)normalitu globalne
library(MVN)
mvn(dec, mvnTest = "royston")

# Jasne, ze to nemoze byt zdruzene normalne (vid tiez Ashton Eaton)
hist(apply(dec, 1, sum), breaks = 100)

# Pozrime si vzajomne vztahy dvojic premennych, ich korelacie a kovariancie
plot(dec, pch = 19, cex = 0.1)
plot(dec[, 1:6], pch = 19, cex = 0.1)
round(cov(dec), 3)
diag(cov(dec))
plot(sqrt(diag(cov(dec))), pch = 19)
round(cor(dec), 3)

# Vypocitajme PCA "nizkourovnovo"
S <- cov(dec)
lam <- eigen(S)$values; lam
U <- eigen(S)$vectors; U

# Zavislost na skalovani premennych - problem, ktory ma PCA v DNA!
# Preskalujme
R <- cor(dec)
lam <- eigen(R)$values; lam
U <- eigen(R)$vectors; U

# Pouzime teraz standardnu fciu prcomp
dec.pca <- prcomp(dec, scale = TRUE)
dec.pca
summary(dec.pca)

# Laktovy diagram
plot(0:10, c(0, summary(dec.pca)$importance[3, ]), type = "b", ylim = c(0, 1))

# Nakuknime dovnutra objektu dec.pca
dec.pca[1:5]

# Ako by sme mohli interpretovat a nazvat prve tri hlavne komponenty?
round(dec.pca$rotation[, 1], 1)
round(dec.pca$rotation[, 2], 1)
round(dec.pca$rotation[, 3], 1)

# Zobrazenie v rovine pomocou PCA, najprv "nizkourovnovo"
plot(dec.pca$x[, 1], dec.pca$x[, 2], type = "n")
text(dec.pca$x[, 1], dec.pca$x[, 2], rownames(dec), cex = 0.7)

# Kedze chapeme o co ide, mozeme pouzit hotovu funkciu
biplot(dec.pca, cex = 0.7)

# Na koniec pochopme o co ide tu: https://en.wikipedia.org/wiki/Eigenface