#setwd("c:/Users/janousova/Documents/Vyuka/Vicerozmerky a AKD/Analyza a klasifikace dat - vyuka/2014-podzim/AKD_predn-05")



######### Priklad - klasifikace #############

X1=matrix(c(2,4,3,12,10,8),3,2) # data pacientu
X1
X2=matrix(c(5,3,4,7,9,5),3,2) # data kontrol
X2
mm1=as.matrix(colMeans(X1,1)) # centroid pacientu - dve m jsou tu pouzity proto, aby bylo patrne, ze je to vektor
mm1
mm2=as.matrix(colMeans(X2,1)) # centroid kontrol - dve m jsou tu pouzity proto, aby bylo patrne, ze je to vektor
mm2


### Fisherova LDA:
S1=var(X1) # kovariancni matice pacientu
S1
S2=var(X2) # kovariancni matice kontrol
S2
#Sw=((nrow(X1)-1)*S1+(nrow(X2)-1)*S2)/(nrow(X1)-1+nrow(X2)-1) # vazeny vypocet Sw je podle knizky Applied Multivariate Statistical Analysis
Sw=S1+S2
Sw
Sw.inv=solve(Sw)
Sw.inv

# vypocet diskriminacniho smeru:
w=Sw.inv%*%(mm1-mm2)
w
w=w*6 # preskalovani na cela cisla
w 
#w=w/(sqrt(sum(w^2))) % normalizace vektoru - pak je hranicni bod primo ve vzdalenosti m od pocatku
#w

m1=t(w)%*%mm1 # projekce centroidu pacientu
m1
# nebo take: m1=t(Sw.inv%*%(mm1-mm2))%*%mm1 *6
m2=t(w)%*%mm2 # projekce centroidu kontrol
m2
m=(m1+m2)/2 # hranicni bod (prah)
# nebo take: m=(1/2)*t(mm1-mm2)%*%Sw.inv%*%(mm1+mm2) *6 # podle (11-25) v knizce Applied Multivariate Statistical Analysis
m

#library(plotrix)
#draw.circle(0,0,m,border="yellow",lty=1,lwd=1) # pokud je pouzit normovany vektor w, je hranicni bod ve vzdalenosti m od pocatku


# vykresleni:
library("MASS")
eqscplot(X1,xlim=c(-3,12),ylim=c(-3,12),col="red",pch=16,cex=1.5) # pacienti (eqscplot je pouzit, aby obe osy mely stejne meritko, aby pak kolmice byly na sebe opravdu kolme)
#plot(X1,xlim=c(-1,12),ylim=c(-1,12),col="red",pch=16,cex=1.5) # pacienti (nestaci pouzit jen stejne rozsahy os, aby mely osy stejne meritko, protoze bile plochy nejsou stejne velke na vsech stranach obrazku)
points(X2,pch=16,cex=1.5) # kontroly
points(t(mm1),pch=3,col="red",cex=1.5,lwd=2) # centroid pacientu
points(t(mm2),pch=3,cex=1.5,lwd=2) # centroid kontrol
x0=as.matrix(c(3.5,9)) # obraz, ktery se ma klasifikovat
points(t(x0),pch=16,col="blue",cex=1.5) # vykresleni obrazu, ktery se ma klasifikovat
lines(c(-3,12),c((-w[1]/w[2])*(-3)+(m/w[2]),(-w[1]/w[2])*12+(m/w[2])),col="blue",lwd=2) # klasifikacni hranice
lines(c(0,0),c(-3,12)) # osa y
lines(c(-3,12),c(0,0)) # osa x

points(0,m/w[2],col="blue") # prusecik hranice s osou y
points(m/w[1],0,col="blue") # prusecik hranice s osou x - nevejde se do obrazku, je tu jen proto, aby bylo videt, jak se pocita

lines(c(mm1[1],mm2[1]),c(mm1[2],mm2[2])) # hranice puli spojnici centroidu
lines(c(-3,12),c((w[2]/w[1])*(-3),(w[2]/w[1])*(12)),col="cyan",lwd=2) # diskriminacni smer

mm=as.matrix(solve(matrix(c(w[1],w[2],w[2],-w[1]),2,2),c(m,0))) # hranicni bod v prostoru (je to prusecik hranice a nove osy - resi se soustava obecnych rovnic techto dvou primek)
mm
t(w)%*%mm # overeni, ze projekci ziskame bod m=13.5
points(t(mm),col="cyan") # vykresleni hranicniho bodu


# klasifikace noveho subjektu:
x0=as.matrix(c(3.5,9))
y0=t(w)%*%x0 # projekce bodu reprezentujiciho novy subjekt
y0
if (y0>=m) {
   group=1
} else {
   group=2
}
group # novy subjekt je klasifikovan jako pacient



######################################################################################################################


# klasifikace - resubstituce:
gt=c(array(1,nrow(X1)),array(2,nrow(X2))) # ground truth (skutecna prislusnost subjektu do skupin)
group=array(0,nrow(X1)+nrow(X2)) # alokace vektoru, kam budou ukladany vysledky klasifikace jednotlivych subjektu
for (i in 1:nrow(X1)) { # klasifikace dat X1 (pacienti)
   if ((t(w)%*%as.matrix(X1[i,]))>=m) {
      group[i]=1
   } else {
      group[i]=2
   }
}
for (i in 1:nrow(X2)) { # klasifikace dat X2 (kontroly)
   if ((t(w)%*%as.matrix(X2[i,]))>=m) {
      group[i+nrow(X1)]=1
   } else {
      group[i+nrow(X1)]=2
   }
}
gt
group
# je patrne, ze je spatne klasifikovany jeden pacient a jeden kontrolni subjekt (je to videt i z obrazku)
# je to ukazka, ze LDA nemusi byt vzdy vyhra - bylo by urcite mozne navrhnout hranici tak, ze by byl spatne klasifikovany jen ten jeden kontrolni subjekt


### odloz-jeden-mimo krizova validace (LOOCV):
gt=c(array(1,nrow(X1)),array(2,nrow(X2))) # ground truth (skutecna prislusnost subjektu do skupin)
group2=array(0,nrow(X1)+nrow(X2)) # alokace vektoru, kam budou ukladany vysledky klasifikace jednotlivych subjektu
# testovaci data budou jednotlivi pacienti:
mm2=as.matrix(colMeans(X2,1)) # centroid kontrol
S2=var(X2) # kovariancni matice kontrol
for (i in 1:nrow(X1)) { 
   x_test=X1[i,] # testovaci subjekt
   X1_train=X1[-i,] # vyber trenovacich dat pacientu bez testovaciho subjektu
   mm1=as.matrix(colMeans(X1_train,1)) # centroid pacientu
   S1=var(X1_train) # kovariancni matice pacientu
   Sw=S1+S2 # suma ctvercu rozptylu uvnitr skupin
   Sw.inv=solve(Sw)
   w=Sw.inv%*%(mm1-mm2) # diskriminacni smer
   m1=t(w)%*%mm1 # projekce centroidu pacientu
   m2=t(w)%*%mm2 # projekce centroidu kontrol
   m=(m1+m2)/2 # hranicni bod (prah)
   y_test=t(w)%*%x_test # projekce bodu reprezentujiciho testovaci subjekt
   if (y_test>=m) {
      group2[i]=1
   } else {
      group2[i]=2
   }
}
# testovaci data budou jednotlive kontrolni subjekty:
mm1=as.matrix(colMeans(X1,1)) # centroid pacientu
S1=var(X1) # kovariancni matice pacientu
for (i in 1:nrow(X1)) { 
   x_test=X2[i,] # testovaci subjekt
   X2_train=X2[-i,] # vyber trenovacich dat kontrol bez testovaciho subjektu
   mm2=as.matrix(colMeans(X2_train,1)) # centroid kontrol
   S2=var(X2_train) # kovariancni matice kontrol
   Sw=S1+S2 # suma ctvercu rozptylu uvnitr skupin
   Sw.inv=solve(Sw)
   w=Sw.inv%*%(mm1-mm2) # diskriminacni smer
   m1=t(w)%*%mm1 # projekce centroidu pacientu
   m2=t(w)%*%mm2 # projekce centroidu kontrol
   m=(m1+m2)/2 # hranicni bod (prah)
   y_test=t(w)%*%x_test # projekce bodu reprezentujiciho testovaci subjekt
   if (y_test>=m) {
      group2[i+nrow(X1)]=1
   } else {
      group2[i+nrow(X1)]=2
   }
}
gt
group2
# spatne klasifikovani dva pacienti a dva kontrolni subjekty