# Ilustracja liniowej analizy dyskryminacyjnej (lda) na krabach i porownanie jej
# z PCA. Kolory sa zgodne z nazwami podpopulacji. Samce oznaczone sa trojkatami.
# LDA policzylem "na piechote" - zakomentarzowane sa wersje standardowe.
# Liniowe Funkcje Dyskryminacyjne = Zmienne Kanoniczne  = kolumny Z
# Kierunki najwiekszego zróżnicowania = wektory kanoniczne =
#                                     = wektory dyskryminacyjne = kolumny  V
# Z = X%*%V
# T = B + W, gdze T - total covariance of data, B - between groups,
# W - within groups. Rozwiazuje uogólnione zadanie własne : B*v = la*W*v

rozkladT=function(clas) {
  n=nrow(X); mX=apply(X,2,mean)
  n0=length(levels(factor(clas)))
  p=ncol(X); W=matrix(0,p,p); B=W
  for(i in 1:n0){
    Xi=X[clas==i,];  ni=nrow(Xi);  mi=apply(Xi,2,mean)
    W= W + (ni-1)*var(Xi);  B= B + ni*(mi-mX)%*%t(mi-mX)
  }
  list(W=W/n,B=B/n)
}

library(MASS); X=as.matrix(crabs[,4:8])
# Normalizowanie danych lub przeksztalcenie log:
#X=t((t(X)-apply(X,2,mean))/apply(X,2,sd)); #X=log(X)
clas=c(rep(1,50),rep(2,50),rep(3,50),rep(4,50))
specie=c(rep("blue",100),rep("orange",100))
sex= c(rep(2,50),rep(1,50),rep(2,50),rep(1,50))
# Na piechotę
lst=rozkladT(clas); W=lst$W; B=lst$B
Uinv=backsolve(chol(W),diag(nrow(W))) # W = t(U)%*%U
E=eigen(t(Uinv)%*%B%*%Uinv);  V=Uinv%*%E$vec; Z=X%*%V

#postscript("lda.eps", onefile=FALSE, horizontal=TRUE, pointsize=10)
par(mfrow=c(1,3),pty="s")
plot(-Z[,1], Z[,2], xlab="LD1", ylab="LD2", col=specie, pch=sex, main="LDA", cex=1.2)

# Inne metody obliczania LDA
#plot(X%*%lda(X,clas)$scal[,1:2], col=specie, pch=sex, main="lda", cex=1.2)

# Ta metoda inaczej skaluje linear discriminant functions
#plot(predict(lda(X,clas))$x, col=specie, pch=sex, main="predict(lda)", cex=1.2)

# Jeszcze inny sposób
#Y=cbind( c(rep(1,50),rep(0,150)),
#         c(rep(0,50),rep(1,50),rep(0,100)),
#         c(rep(0,100),rep(1,50),rep(0,50)))#,rep(0,200) )
# Inaczej zrobiony Y: library(nnet); Y=class.ind(cl)
#CC=cancor(X,Y); Z.ca=X%*%as.matrix(CC$xcoe)
#plot(Z.ca[,1:2], xlab="ZmKanon1",ylab="ZmKanon2", col=clas, main="zm kanoniczne")

# PCA dla porównania
plot(prcomp(X)$x[,1:2], col=specie, pch=sex, main="PCA", cex=1.2)
plot(prcomp(X,scale=TRUE)$x[,1:2], col=specie, pch=sex, main="PCA scaled", cex=1.2)
#graphics.off()
print("Ranking udzialu pierwotnych cech w zmiennych kanonicznych")
print(cor(Z[,1:3],X)^2); cat("\n")
print("Wektory dyskryminacyjne * sqrt(diag(W))")
A=t(V[,1:3]*sqrt(diag(W))); colnames(A)=colnames(X)
print(A)