##=============================
# Kovariancna, korelacna matica
#------------------------------

#install.packages("corrplot") # toto spustit iba pri prvom pouziti
library(corrplot)   # kvoli funkcii corrplot

# Otcovia a synovia (Pearson, cca. 1900)
#------------------

# povodne data, v palcoch
Z <- read.table("father.son.csv")

# data v centimetroch
X <- Z*2.54
colnames(X) <- c("otcovia", "synovia")

# prvych par udajov
head(X)

plot(X)

# Kovariancna matica - rozptylenosti, zavislosti
cov(X)

# kontrola manualneho vzorca
n <- nrow(X)
P <- diag(n) - matrix(1, n, n)/n
S <- t(as.matrix(X)) %*% P %*% as.matrix(X) / (n-1)
S
S - cov(X)   # to iste

# Korelacna matica - korelacie (zavislosti)
cor(X)

# kontrola manualneho vzorca
s <- 1/sqrt(diag(S))
R <- diag(s) %*% S %*% diag(s)
R
cor(X) - R   # to iste

# vykreslenie kor. matice
corrplot(cor(X))

# Wine tasting
#-------------

# udaje o zlozeniach vin

M <- read.table("combwines.csv", sep=',', header=TRUE)
X.wine <- M[, -c(12,13)]
lab <- M[,13]   # kvalita
lac <- M[,12]   # biele/cervene

# kovariancna matica - rozptylenosti, zavislosti
S.wine <- cov(X.wine)
S.wine

# korelacna matica - korelacie (zavislosti)
R.wine <- cor(X.wine)
R.wine

# vykreslenie kor. matice
corrplot(R.wine)

##================================
# Viacrozmerne normalne rozdelenie
#---------------------------------

# N(0, Sigma)
k <- 2
n <- 1000
A <- diag(k) - matrix(1, nrow=k, ncol=k)/sqrt(k)
Sigma <- A %*% t(A)

# install.packages("rgl")
# install.packages("mvtnorm")
library(rgl)
library(mvtnorm)

# vykreslenie hustoty
x <- seq(-2, 2, by=0.1)
y <- seq(-2, 2, by=0.1)
z <- matrix(nrow=length(x), ncol=length(y))
for(i in 1:length(x)) {
   for(j in 1:length(y)) {
      z[i,j] <- dmvnorm(c(x[i], y[j]), sigma=Sigma)
   }
}

open3d()
persp3d(x, y, z, col="red")

# vrstevnice
library(ellipse)
plot(ellipse(Sigma, t=5), type="l")
for(i in 1:4) {
   lines(ellipse(Sigma, t=i))
}

# spektralny rozklad
res <- eigen(Sigma)
lam <- res$values
print(lam)
Q <- res$vectors
print(Q)

# vlastne vektory Sigma
Q2 <- Q %*% diag(sqrt(lam)) # aby dlzky zodpovedali sqrt(lam)
lines(c(0, 3*Q2[1,1]), c(0, 3*Q2[2,1]), col="red", lwd=2)
lines(c(0, 3*Q2[1,2]), c(0, 3*Q2[2,2]), col="red", lwd=2)

# nahodny vyber z N(0, Sigma)
xnorm <- rmvnorm(200, sigma=Sigma)
points(xnorm, pch=16, cex=0.5)

##===================
# PCA
#--------------------

# Otcovia a synovia

# centrovanie
mX <- apply(X, 2, mean)
X <- X - matrix(mX, nrow=nrow(X), ncol=ncol(X), byrow=TRUE)
apply(X, 2, mean)
lims <- c(min(X), max(X))
plot(X, xlim=lims, ylim=lims)

# spektralny rozklad S a vykreslenie
resS <- eigen(S)
lam <- resS$values
Q <- resS$vectors
Q2 <- Q %*% diag(sqrt(lam)) # aby dlzky zodpovedali sqrt(lam)

# vlastne vektory
lines(c(0, 3*Q2[1,1]), c(0, 3*Q2[2,1]), col="red", lwd=3)
lines(c(0, 3*Q2[1,2]), c(0, 3*Q2[2,2]), col="red", lwd=3)
   # najviac rozptylene: v smere q1

# Wine tasting
#-------------

# vykresleny nahodny vyber (prilis vela bodiek je inak)
n <- nrow(X.wine)
rndrows <- sample(n, 1000)
Xrnd <- X.wine[rndrows,]
plot(Xrnd, col="darkgrey", cex=0.5)
plot(Xrnd, col=c("#227700", "#cc0000")[lab][rndrows], cex=0.5)

# PCA
wine.pca <- prcomp(X.wine, scale=TRUE)
Y <- wine.pca$x[,c(1,2)]
plot(Y, cex=0.5, pch=16, col="darkgrey")
plot(Y, cex=0.5, pch=16, col=c("#227700", "#cc0000")[lab])

# dokreslenie kvalitnych vin
plot(Y, cex=0.5, pch=16, col="#999999")
points(Y[lac>=7, ], cex=0.5, pch=16, col="#000099")

# Manualne:

# naskalovanie
mX.wine <- apply(X.wine, 2, mean)
X.wine <- X.wine - matrix(mX.wine, nrow=nrow(X.wine), ncol=ncol(X.wine), byrow=TRUE)

s.wine <- 1/sqrt(diag(S.wine))
X.wine <- as.matrix(X.wine) %*% diag(s.wine)

S.wine <- cov(X.wine)   # S pre naskalovane data

# spektralny rozklad
res.wine <- eigen(S.wine)
lam.wine <- res.wine$values
Q.wine <- res.wine$vectors

# PCA:
Y.m <- as.matrix(X.wine) %*% Q.wine
plot(Y.m[, 1:2])

par(mfrow=c(1,2))
plot(Y, cex=0.5, pch=16, col=c("#227700", "#cc0000")[lab], main="Rko")
plot(Y.m, cex=0.5, pch=16, col=c("#227700", "#cc0000")[lab], main="manualne")
dev.off()

# namiesto q2 si Rko zvolilo -q2:
Q.wine[, 2] <- - Q.wine[, 2]
Y.m <- as.matrix(X.wine) %*% Q.wine
   # takto zmenit/nezmenit prvy/druhy stlpec
   # ... podla toho, ktora suradnica vysla "naopak"

par(mfrow=c(1,2))
plot(Y, cex=0.5, pch=16, col=c("#227700", "#cc0000")[lab], main="Rko")
plot(Y.m, cex=0.5, pch=16, col=c("#227700", "#cc0000")[lab], main="manualne")
dev.off()

#==================
# Linearna regresia
#------------------

# Otcovia a synovia
#------------------

Z <- read.table("father.son.csv")
data <- Z*2.54
colnames(data) <- c("otcovia", "synovia")

# vykreslenie
plot(data, xlab="vyska otca", ylab="vyska syna",
     cex.axis=1.25, cex.lab=1.5, pch=16, cex=0.5)

# regresna priamka
m <- lm(synovia ~ otcovia, data=data)
abline(m$coef, lwd=3, col="blue")
m
   # ak je otec vyssi o 1 cm, tak syn je v priemere vyssi o 0.5 cm

# manualne
y <- data$synovia
X <- cbind(1, data$otcovia)
head(X)

beta <- solve(t(X) %*% X) %*% t(X) %*% y
beta
m$coef   # to iste

# Boston housing
#---------------

#install.packages("MASS") # toto spustit iba pri prvom pouziti
library(MASS)

data(Boston)
# hodnoty domov v primestskych castiach Bostonu (v 1000 USD)
# medv ... medianova hodnota domu
# crim ... miera kriminality (per capita)
# nox ... koncetracia oxidov dusika
# rm ... priemerny pocet izieb
# dis ... priemerna vzdialenost od oblasti s velkou koncentraciu zamestnancov
# tax ... miera zdanenia nehnutelnosti (na 10 000 USD)

head(Boston)

m2 <- lm(medv ~ crim + nox + rm + dis + tax, data=Boston)
m2
   # zvysenie miery kriminality o 1 sposobuje odhadovany pokles medianovej hodnoty 
   # o 0.15, t.j. o 150 USD
   # apod.