# Predmet: "Analyza zhlukov a klasifikacia dat"
# Studijny program: "Pravdepodobnost a matematicka statistika"
# Vyucujuci: Radoslav Harman, KAMS FMFI UK Bratislava
#
# Linearna diskriminacna analyza
# Data Tibet a ciastocne aj analyza je zalozena na kapitole 7 v knihe
# Everitt, B: R and S+ Companion to Multivariate Analysis

Tibet <- read.table("http://www.iam.fmph.uniba.sk/ospm/Harman/data/tibet.txt", header = TRUE)
dim(Tibet); n <- nrow(Tibet)
print(Tibet)

# Data zodpovedaju piatim meraniam (vysvetlujucim premennym, prediktorom, crtam)
# dvoch typov (kategorii, tried) lebiek (objektov):
# 17 lebiek je typu I (tibetsky kmen Sikkim) a 15 je typu II (tibetsky kmen Khan). 
# Vysvetlujuce premenne: Length, Breadth, Height, Face height, Face breadth

attach(Tibet)

# Zobrazme si data v klasickych 2D priemetoch aj priemetoch urcenych cez PCA
plot(Tibet[, -6], col = Type, pch = 19, asp = 1)
plot(as.data.frame(prcomp(Tibet[, -6])$x), col = Type, pch = 19,
     xlim = c(-30, 30), ylim = c(-30, 30), asp = 1)
print(prcomp(Tibet[, -6]))

# Extrahujme z treningovych dat len crty (vyhodime spravnu klasifikaciu)
x1 <- Tibet[Type == 1, -6]
x2 <- Tibet[Type == 2, -6]

# Vyskusajme, ci je akceptovatelna zdruzena normalita vysvetlujucich premennych
# pre kazdu skupinu (toto nie je nutna podmienka, aby bol LDA/QDA klasifikator uzitocny)
library(MVN)
mvn(x1, mvnTest = "royston")
mvn(x2, mvnTest = "royston")

# Vyskusajme tiez ci je akceptovatelna hypoteza o rovnosti kovariancnych matic
# (nie nutna podmienka na to, aby bol LDA klasifikator uzitocny)
library(MVTests)
BoxM(Tibet[, -6], Type)

# Vypocitame obdhad strednych hodnot vektorov crt pre oba kmene zvlast
m1 <- apply(x1, 2, mean); print(m1)
m2 <- apply(x2, 2, mean); print(m2)

# Vypocitajme poolovanu vyberovu kovariancu maticu
n1 <- dim(x1)[1]; n2 <- dim(x2)[1]
S123 <- ((n1 - 1) * var(x1) + (n2 - 1) * var(x2)) / (n1 + n2 - 2)

# Vykonajme test rovnosti strednych hodnot z predmetu VSA minuly semester
T2 <- (n1 * n2) / (n1 + n2) * t(m1 - m2) %*% solve(S123) %*% (m1 - m2)
Fstat <- (n1 + n2 - 5 - 1) * T2 / ((n1 + n2 - 2) * 5) # Tu ma Everitt drobnu chybicku
pvalue <- 1 - pf(Fstat, 5, 26)
print(pvalue)
# Vidime, ze stredne hodnoty su statisticky vyznamne posunute

# Vypocet linearnej diskriminacnej funkcie (formalne, l(x)=a^T*x)
a <- solve(S123) %*% (m1 - m2)
print(a)

z12 <- as.numeric((m1 %*% a + m2 %*% a) / 2)
print(z12)

# To je to iste ako vzorec z prednasky pre b, aha:
b <- 0.5 * (t(m2) %*% solve(S123) %*% m2 - t(m1) %*% solve(S123) %*% m1)
print(b)

# Teraz si ukazeme ako toto urobime rychlejsie s uz implementovanou funkciou v Rku 
library(MASS); help(lda)
dis <- lda(Type ~ Length + Breadth + Height + Fheight + Fbreadth,
           data = Tibet, prior = c(1/2, 1/2))
# parametere prior su presne q1 a q2 z prednasky. Berieme L(2|1) = L(1|2)
# Ine straty L(2|1) a L(1|2) sa daju do procedury vlozit umelou zmenou q1 a q2 

print(dis$means)
print(dis$scaling)
# Toto sa zda byt ine ako vyslo nam z teorie, avsak diskriminaciu to robi taku istu:

print(dis$scaling / a)
print((t(dis$scaling) %*% (m1 + m2) / 2)/z12)

# Urobme odhady presnosti klasifikatora

# Resubstitucia (in-sample method) a zobrazenie
# Nizkourovnovo
cls <- 1.5 - sign(as.matrix(Tibet[, -6]) %*% a - z12) / 2
print(cls)

# Pouzitim predict.lda
Tibet.pred <- predict(dis, Tibet[, -6], prior = c(1/2, 1/2), method = "plug-in")
print(Tibet.pred)

# Pozrime sa, ci je klasifikcia naozaj rovnaka
cbind(cls, Tibet.pred$class)

# Zobrazme si omyly v rovine prvych dvoch PCs
x1.pca <- prcomp(Tibet[, -6])$x[, 1]
x2.pca <- prcomp(Tibet[, -6])$x[, 2]

plot(x1.pca, x2.pca, pch = 19, 
     col = as.numeric(Type), cex = 1 + 2*abs(cls - Tibet[, 6]))
# Resubsitucne odhady chybnych klasifikacii su teda 3/17, 3/15

# Skusime aj odhad pomocou Mahalanobisovej vzdialenosti, ktory je specificky pre LDA
alpha2 <- t(m1 - m2) %*% solve(S123) %*% (m1 - m2)
k <- 1 # Zamerne, aby bol jasny suvis so vzorcami z prednasky
print(1 - pnorm((log(k) + alpha2 / 2) / sqrt(alpha2))) # Odhad P(1|2)
print(pnorm((log(k) - alpha2 / 2) / sqrt(alpha2))) # Odhad P(2|1)
# Odhady P(1|2) a P(2|1) pomocou M-vzdialenosti su tu rovnake,
# Vsimnite si, ze podla vzorcov je to tak vzdy ked k=1

# Urobme aj krizovu validaciu typu "leave one out"
conf.mat <- matrix(0, ncol = 2, nrow = 2) # Toto bude nenormalizovana confusion matrix

for (i in 1:nrow(Tibet)) {
  
  dis <- lda(Type ~ Length + Breadth + Height + Fheight + Fbreadth,
             data = Tibet[-i, ], prior = c(1/2, 1/2))
  m1 <- as.vector(dis$means[1, ]); m2 <- as.vector(dis$means[2, ])
  a <- as.vector(dis$scaling); b <- as.numeric(t(dis$scaling) %*% (m1 + m2) / 2)
  
  # Vypocitame ako LDA natrenovana na datach Tibet[-i, ]
  # bude klasifikovat vynechany objekt Tibet[i, ] 
  cls <- 1.5 - sign(b - sum(as.vector(Tibet[i, 1:5]) * a)) / 2
  conf.mat[Tibet[i, 6], cls] <- conf.mat[Tibet[i, 6], cls] + 1
  if (Tibet[i, 6] != cls) text(x1.pca[i], x2.pca[i], labels = "X", col = "green", cex = 2)
}

confusion.info <- function(nnC) {
  # Nenormalizovana aj normalizovana confusion matrix a odhady roznych charakteristik
  # nnC ... nenormalizovana kontingencna tabulka klasifikacie (confusion matrix)
  print("Non-normalized confusion matrix"); print(nnC)
  C <- nnC/sum(nnC); print("Normalized confusion matrix"); print(round(C, 4))
  print(c("Sensitivity", round(C[1, 1]/(C[1, 1] + C[1, 2]), 4)))
  print(c("Specificity", round(C[2, 2]/(C[2, 1] + C[2, 2]), 4)))
  print(c("Accuracy", round(C[2, 2] + C[1, 1], 4)))
  PY <- (C[1, 1]*C[2, 2] - C[1, 2]*C[2, 1]) /
    sqrt(sum(C[1, ])*sum(C[2, ])*sum(C[, 1])*sum(C[, 2]))
  print(c("Pearson-Yule", round(PY, 4)))
}

confusion.info(conf.mat)

# Nakoniec si este mozeme nakreslit resubstitucny odhad ROC krivky

# Radsej si opat prepocitajme vektor a urcujuci diskrimonacnu funkciu:
a <- lda(Type ~ Length + Breadth + Height + Fheight + Fbreadth,
           data = Tibet, prior = c(1/2, 1/2))$scaling

thr <- seq(13, 20, by = 0.01); n.thr <- length(thr)
senz <- spec <- rep(0, n.thr) 
for (i.thr in 1:n.thr) {

  # Triedy, do ktorych klasifikator zaradi objekty pre prah thr[i.thr]
  cls <- 1.5 - sign(thr[i.thr] - as.matrix(Tibet[, -6]) %*% a) / 2
  
  # Aktualna confusion matrix
  C <- matrix(0, ncol = 2, nrow = 2)
  for (i in 1:n) {
    C[Type[i], cls[i]] <- C[Type[i], cls[i]] + 1
  }
  
  senz[i.thr] <- C[1, 1]/(C[1, 1] + C[1, 2])
  spec[i.thr] <- C[2, 2]/(C[2, 1] + C[2, 2])
}

plot(1 - spec, senz, type = "l", xlim = c(0, 1), ylim = c(0, 1), lwd = 2); grid()