MyLDA <- function(X,C) # X: Datenmatrix (n Beobachtungen = n Zeilen) # C: Vektor von Klassenzugehörigkeiten # in Form von Zahlen in {1,2,...,L} fuer ein L, # wobei alle L Werte vorkommen! { # Dimensionen: tmp <- dim(X) n <- tmp[1] p <- tmp[2] # Anzahl Gruppen und deren Groessen: N <- table(C) L <- max(C) # Mittelwert aller Trainingsvektoren: Xbar <- apply(X,2,mean) # Totale Streuung: Xc <- X - matrix(rep(1,n),ncol=1) %*% matrix(Xbar,nrow=1) SSP.total <- t(Xc) %*% Xc # Gruppenweise Mittelwerte: Mhat <- matrix(0,nrow=L,ncol=p) for (y in 1:L) { Mhat[y,] <- apply(X[C==y,],2,mean) } Mhat # Streuung innerhalb der Gruppen: str(X) str(Mhat) X.intra <- X - Mhat[C,] SSP.intra <- t(X.intra) %*% X.intra SSP.intra # Streuung zwischen den Gruppen: tmp <- Mhat - matrix(rep(Xbar,times=L),nrow=L,byrow=TRUE) SSP.inter <- t(tmp) %*% diag(N) %*% tmp SSP.inter # Berechnung von "Sigma-hut hoch -1/2": tmp <- eigen(SSP.intra/(n-L)) Shm12 <- tmp$vectors %*% diag(1/sqrt(tmp$values)) %*% t(tmp$vectors) # Streuung zwischen den Gruppen fuer die mit # Shm12 transformierten Trainingsvektoren: Gamma <- Shm12 %*% SSP.inter %*% Shm12 tmp <- eigen(Gamma) # Matrix der Diskriminanzvektoren (Spalten): B <- Shm12 %*% tmp$vectors # Transformierte Datenmatrix: X.trafo <- X %*% B # Graphische Darstellung fuer d=2: par(mai=c(0.5,0.5,0.1,0.1)) plot(X.trafo[,1], X.trafo[,2],pch=19,xlab="",ylab="") Farben <- c("red","green","blue","magenta","cyan","yellow") if (L > 7) { Farben <- rep(Farben,times=ceiling(L/7)) } Muster <- c(22,23,21,24,25) if (L > 5) { Muster <- rep(Muster,times=ceiling(L/5)) } for (y in (1:L)) { points(X.trafo[C==y,1],X.trafo[C==y,2], pch=Muster[y],col="black",bg=Farben[y]) } return(list(L=L,N=N,Xbar=Xbar,Mhat=Mhat, SSP.total=SSP.total, SSP.intra=SSP.intra, SSP.inter=SSP.inter, B=B)) }