sad_11_lda

2276 days ago by macieksk

# Przyklady Piotr Pokarowski # http://www.mimuw.edu.pl/~pokar/StatystykaII/ 
       
# http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/lda.R # Opis metody LDA # http://miner.chem.purdue.edu/Lectures/Lecture10.pdf 
       
# 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 zroznicowania = 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 uogolnione zadanie wlasne : 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)) 
       
head(crabs) head(X) head(clas) 
       
  sp sex index   FL  RW   CL   CW  BD
1  B   M     1  8.1 6.7 16.1 19.0 7.0
2  B   M     2  8.8 7.7 18.1 20.8 7.4
3  B   M     3  9.2 7.8 19.0 22.4 7.7
4  B   M     4  9.6 7.9 20.1 23.1 8.2
5  B   M     5  9.8 8.0 20.3 23.0 8.2
6  B   M     6 10.8 9.0 23.0 26.5 9.8
    FL  RW   CL   CW  BD
1  8.1 6.7 16.1 19.0 7.0
2  8.8 7.7 18.1 20.8 7.4
3  9.2 7.8 19.0 22.4 7.7
4  9.6 7.9 20.1 23.1 8.2
5  9.8 8.0 20.3 23.0 8.2
6 10.8 9.0 23.0 26.5 9.8
[1] 1 1 1 1 1 1
# Na piechote 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) plot(-Z[,1], Z[,2], xlab="LD1", ylab="LD2", col=specie, pch=sex, main="LDA", cex=1.2) dev.off() 
       
null device 
          1 
png(w=1000) layout(matrix(1:2,1)) # 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) dev.off() 
       
null device 
          1 
# Jeszcze inny sposob 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") dev.off() 
       
null device 
          1 
# PCA dla porownania 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) dev.off() 
       
null device 
          1 

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) 
       
[1] "Ranking udzialu pierwotnych cech w zmiennych kanonicznych"
             FL         RW         CL          CW         BD
[1,] 0.19189666 0.14725468 0.07250842 0.040006096 0.17226895
[2,] 0.01070583 0.08843377 0.01969958 0.007358558 0.02625701
[3,] 0.36228633 0.22782251 0.32430869 0.354517454 0.24689353

[1] "Wektory dyskryminacyjne * sqrt(diag(W))"

             FL        RW        CL         CW        BD
[1,]  4.8134222  1.422677  1.259730 -11.473183  4.167556
[2,] -0.6044627 -3.505707  7.357538  -4.871444  1.584833
[3,]  5.1615777 -1.038115 -4.572525   4.957266 -3.954936
###################################333 ############# QDA - Quadratic_discriminant_analysis 
       
# http://www.mimuw.edu.pl/~pokar/StatystykaII/PREDYKCJA/CrossValKlasCrabs.R library(MASS); library(nnet) X=as.matrix(crabs[,4:8]); y=c(rep(1,50),rep(2,50),rep(3,50),rep(4,50)) pe1=sample(50); pe2=sample(50)+50; pe3=sample(50)+100; pe4=sample(50)+150 ppp.qda=rep(NA,5); ppp.nnet=rep(NA,5) head(crabs) head(X) 
       
for(i in 1:5){ ind=((i-1)*10+1):(i*10) te1=pe1[ind]; te2=pe2[ind]; te3=pe3[ind]; te4=pe4[ind] te=c(te1,te2,te3,te4); tr=(1:200)[-te] K=qda(X[tr,],y[tr]) pred=predict(K,X[te,])$class ppp.qda[i] = mean( y[te] == pred ) S=nnet(X[tr,],class.ind(y[tr]),maxit=200, trace=F,size=16,softmax=T) pred2=apply(predict(S,X[te,]),1,which.max) ppp.nnet[i] = mean( y[te] == pred2 ) } cat("\nP-stwo poprawnej predykcji qda:\n",round(ppp.qda,2),"\n") cat("P-stwo poprawnej predykcji nnet:\n",round(ppp.nnet,2),"\n") 
       
P-stwo poprawnej predykcji qda:
 0.88 0.98 0.98 0.95 0.88 
P-stwo poprawnej predykcji nnet:
 0.88 0.9 0.98 0.92 0.98 
# Zrobcie krzywe ROC dla powyzszych modeli # tak jak w #http://sage.mimuw.edu.pl/home/pub/61 -- sadII_lab8_japonicus_genome