################################################### ### chunk number 1: banknotePlot ################################################### options(width = 68) library(dr) library(RColorBrewer) data(banknote) mypalette <- brewer.pal(3, "Set1")[c(3,1)] klassen <- factor(ifelse(banknote[,7], "Falsch", "Echt")) pairs(banknote[,-7], col = mypalette[klassen], main = "banknote-Datensatz", pch = 20) ################################################### ### chunk number 2: banknoteAusw ################################################### banknote.ausw <- banknote[,c("Bottom", "Top")] plot(banknote.ausw, pch = 20, col = mypalette[klassen], main = "banknote-Datensatz - ausgewählte Variablen") ################################################### ### chunk number 3: LDA ################################################### mu1 <- colMeans(banknote.ausw[klassen == "Echt",]) mu2 <- colMeans(banknote.ausw[klassen == "Falsch",]) Sigma.inv <- solve(cov(banknote.ausw)) fun <- function(x1, mu1, mu2, Sigma.inv, N1, N2){ f <- function(x2, x1, mu1, mu2, Sigma.inv, N1, N2){ (t(c(x1, x2)) %*% Sigma.inv %*% (mu2 - mu1) - 0.5*(t(mu2) %*% Sigma.inv %*% mu2 - t(mu1) %*% Sigma.inv %*% mu1) - log(N1/N2)) } uniroot(f, c(0, 15), x1 = x1, mu1 = mu1, mu2 = mu2, Sigma.inv = Sigma.inv, N1 = N1, N2 = N2)$root } bottom.vec.lda <- c(7, 13) top.vec.lda <- sapply(bottom.vec.lda, fun, mu1 = mu1, mu2 = mu2, Sigma.inv = Sigma.inv, N1 = 100, N2 = 100) ################################################### ### chunk number 4: LDAplot ################################################### plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranke für LDA") lines(bottom.vec.lda, top.vec.lda, lwd = 2, lty = 1, col = mypalette[1]) lines(bottom.vec.lda, top.vec.lda, lwd = 2, lty = 2, col = mypalette[2]) text(8, 8, "Echt", font = 2, col = mypalette[1], cex = 2) text(12, 12, "Falsch", font = 2, col = mypalette[2], cex = 2) ################################################### ### chunk number 5: DLDA ################################################### Sigma.inv <- matrix(0, ncol = 2, nrow = 2) Sigma.inv[1,1] <- 1/var(banknote.ausw[,"Bottom"]) Sigma.inv[2,2] <- 1/var(banknote.ausw[,"Top"]) fun <- function(x1, mu1, mu2, Sigma.inv, N1, N2){ f <- function(x2, x1, mu1, mu2, Sigma.inv, N1, N2){ (t(c(x1, x2)) %*% Sigma.inv %*% (mu2 - mu1) - 0.5*(t(mu2) %*% Sigma.inv %*% mu2 - t(mu1) %*% Sigma.inv %*% mu1) - log(N1/N2)) } uniroot(f, c(0, 15), x1 = x1, mu1 = mu1, mu2 = mu2, Sigma.inv = Sigma.inv, N1 = N1, N2 = N2)$root } bottom.vec.dlda <- c(7, 13) top.vec.dlda <- sapply(bottom.vec.dlda, fun, mu1 = mu1, mu2 = mu2, Sigma.inv = Sigma.inv, N1 = 100, N2 = 100) ################################################### ### chunk number 6: DLDAplot ################################################### plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranke für DLDA") lines(bottom.vec.dlda, top.vec.dlda, lwd = 2, lty = 1, col = mypalette[1]) lines(bottom.vec.dlda, top.vec.dlda, lwd = 2, lty = 2, col = mypalette[2]) text(8, 8, "Echt", font = 2, col = mypalette[1], cex = 2) text(12, 12, "Falsch", font = 2, col = mypalette[2], cex = 2) ################################################### ### chunk number 7: QDA ################################################### Sigma1 <- cov(banknote.ausw[klassen == "Echt",]) Sigma1.inv <- solve(Sigma1) Sigma2 <- cov(banknote.ausw[klassen == "Falsch",]) Sigma2.inv <- solve(Sigma2) fun <- function(x1, mu1, mu2, Sigma1, Sigma1.inv, Sigma2, Sigma2.inv, N1, N2, int){ f <- function(x2, x1, mu1, mu2, Sigma1, Sigma1.inv, Sigma2, Sigma2.inv, N1, N2){ x <- c(x1, x2) D1 <- (- 0.5*log(det(Sigma1)) - 0.5*t(x - mu1) %*% Sigma1.inv %*% (x - mu1) + log(N1/(N1 + N2))) D2 <- (- 0.5*log(det(Sigma2)) - 0.5*t(x - mu2) %*% Sigma2.inv %*% (x - mu2) + log(N2/(N1 + N2))) return(D2 - D1) } uniroot(f, int, x1 = x1, mu1 = mu1, mu2 = mu2, Sigma1 = Sigma1, Sigma1.inv = Sigma1.inv, Sigma2 = Sigma2, Sigma2.inv = Sigma2.inv, N1 = 100, N2 = 100)$root } bottom.vec <- seq(7, 13, 0.1) top.vec <- sapply(bottom.vec, fun, mu1 = mu1, mu2 = mu2, Sigma1 = Sigma1, Sigma1.inv = Sigma1.inv, Sigma2 = Sigma2, Sigma2.inv = Sigma2.inv, N1 = 100, N2 = 100, int = c(5, 15)) top.fun.qda <- splinefun(bottom.vec, top.vec, method = "natural") bottom.vec.qda <- seq(7, 13, 0.01) ################################################### ### chunk number 8: QDAplot ################################################### plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranke für QDA") lines(bottom.vec.qda, top.fun.qda(bottom.vec.qda), lwd = 2, lty = 1, col = mypalette[1]) lines(bottom.vec.qda, top.fun.qda(bottom.vec.qda), lwd = 2, lty = 2, col = mypalette[2]) text(8, 8, "Echt", font = 2, col = mypalette[1], cex = 2) text(12, 12, "Falsch", font = 2, col = mypalette[2], cex = 2) ################################################### ### chunk number 9: DQDA ################################################### Sigma1 <- diag(apply(banknote.ausw[klassen == "Echt",], 2, var)) Sigma1.inv <- solve(Sigma1) Sigma2 <- diag(apply(banknote.ausw[klassen == "Falsch",], 2, var)) Sigma2.inv <- solve(Sigma2) fun <- function(x1, mu1, mu2, Sigma1, Sigma1.inv, Sigma2, Sigma2.inv, N1, N2, int){ f <- function(x2, x1, mu1, mu2, Sigma1, Sigma1.inv, Sigma2, Sigma2.inv, N1, N2){ x <- c(x1, x2) D1 <- (- 0.5*log(det(Sigma1)) - 0.5*t(x - mu1) %*% Sigma1.inv %*% (x - mu1) + log(N1/(N1 + N2))) D2 <- (- 0.5*log(det(Sigma2)) - 0.5*t(x - mu2) %*% Sigma2.inv %*% (x - mu2) + log(N2/(N1 + N2))) return(D2 - D1) } uniroot(f, int, x1 = x1, mu1 = mu1, mu2 = mu2, Sigma1 = Sigma1, Sigma1.inv = Sigma1.inv, Sigma2 = Sigma2, Sigma2.inv = Sigma2.inv, N1 = N1, N2 = N2)$root } bottom.vec <- seq(7, 13, 0.1) top.vec <- sapply(bottom.vec, fun, mu1 = mu1, mu2 = mu2, Sigma1 = Sigma1, Sigma1.inv = Sigma1.inv, Sigma2 = Sigma2, Sigma2.inv = Sigma2.inv, N1 = 100, N2 = 100, int = c(0, 20)) top.fun.dqda <- splinefun(bottom.vec, top.vec, method = "natural") bottom.vec.dqda <- seq(7, 13, 0.01) ################################################### ### chunk number 10: DQDAplot ################################################### plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranke für DQDA") lines(bottom.vec.dqda, top.fun.dqda(bottom.vec.dqda), lwd = 2, lty = 1, col = mypalette[1]) lines(bottom.vec.dqda, top.fun.dqda(bottom.vec.dqda), lwd = 2, lty = 2, col = mypalette[2]) text(8, 8, "Echt", font = 2, col = mypalette[1], cex = 2) text(12, 12, "Falsch", font = 2, col = mypalette[2], cex = 2) ################################################### ### chunk number 11: DAvgl ################################################### mypalette1 <- brewer.pal(4, "Dark2") plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranken") lines(bottom.vec.lda, top.vec.lda, lwd = 2, lty = 1, col = mypalette1[1]) lines(bottom.vec.dlda, top.vec.dlda, lwd = 2, lty = 1, col = mypalette1[2]) lines(bottom.vec.qda, top.fun.qda(bottom.vec.qda), lwd = 2, lty = 1, col = mypalette1[3]) lines(bottom.vec.dqda, top.fun.dqda(bottom.vec.dqda), lwd = 2, lty = 1, col = mypalette1[4]) legend("topright", legend = c("LDA", "DLDA", "QDA", "DQDA"), fill = mypalette1) ################################################### ### chunk number 12: GPLSplot ################################################### library(gpls) fit <- gpls(as.integer(klassen)-1 ~ Bottom + Top, data = banknote.ausw) bottom.vec.gpls <- c(7, 13) top.vec.gpls <- (0.5 - fit$coeff[1] - fit$coeff[2]*bottom.vec.gpls)/fit$coeff[3] plot(banknote.ausw, col = mypalette[klassen], lwd = 2, main = "Diskriminanzschranke für GPLS") lines(bottom.vec.gpls, top.vec.gpls, lwd = 2, lty = 1, col = mypalette[1]) lines(bottom.vec.gpls, top.vec.gpls, lwd = 2, lty = 2, col = mypalette[2]) text(8, 8, "Echt", font = 2, col = mypalette[1], cex = 2) text(12, 12, "Falsch", font = 2, col = mypalette[2], cex = 2) ################################################### ### chunk number 13: kNNplot ################################################### library(class) fit <- knn(banknote.ausw, banknote.ausw, cl = klassen, k = 5) plot(banknote.ausw, col = mypalette[klassen], type = "n", lwd = 2, main = "Klassifikation mittels 5-NN") text(banknote.ausw, as.character(fit), cex = 0.7, col = mypalette[klassen]) ################################################### ### chunk number 14: tunekNN ################################################### library(e1071) knn.obj <- tune.knn(banknote.ausw, klassen, k = 1:10, tunecontrol = tune.control(sampling = "boot")) summary(knn.obj) ################################################### ### chunk number 15: tunekNNplot ################################################### plot(knn.obj, main = "Optimales k für banknote-Daten") ################################################### ### chunk number 16: NNplot ################################################### library(nnet) fit.nn <- nnet(klassen ~ Bottom + Top, data = banknote.ausw, size = 2) pred.nn <- predict(fit.nn, banknote.ausw, type = "class") plot(banknote.ausw, col = mypalette[klassen], type = "n", lwd = 2, main = "Klassifikation mittels NN") text(banknote.ausw, pred.nn, cex = 0.7, col = mypalette[klassen]) ################################################### ### chunk number 17: RFplot ################################################### library(rpart) fit <- rpart(klassen ~ Bottom + Top, data = banknote.ausw, method = "class", control = rpart.control(minsplit = 2)) plot(fit, uniform = TRUE, margin = 0.1, branch = 1, main = "Entscheidungsbaum für Echt vs. Falsch") text(fit, use.n = TRUE, all = TRUE, fancy = TRUE, fwidth = 0.7, fheight = 0.9) ################################################### ### chunk number 18: SVMplot ################################################### library(e1071) data.svm <- data.frame(klassen = klassen, banknote.ausw) fit <- svm(klassen ~ Top + Bottom, data = data.svm) plot(fit, data.svm, symbolPalette = mypalette[c(1,2)], svSymbol = "o", grid = 200, col = brewer.pal(3, "Pastel1")[c(3,1)]) ################################################### ### chunk number 19: SVMplot1 ################################################### fit1 <- svm(klassen ~ Top + Bottom, data = data.svm, gamma = 5) plot(fit1, data.svm, symbolPalette = mypalette[c(1,2)], svSymbol = "o", grid = 200, col = brewer.pal(3, "Pastel1")[c(3,1)]) ################################################### ### chunk number 20: tuneSVM ################################################### svm.obj <- tune.svm(banknote.ausw, klassen, gamma = c(0.2, 0.5, 1, 2, 5), cost = 2^(1:4)) summary(svm.obj) ################################################### ### chunk number 21: tuneSVMplot ################################################### plot(svm.obj, main = "Optimale SVM-Parameter radial-Kern") ################################################### ### chunk number 22: SVMplot2 ################################################### set.seed(13) fit2 <- svm(klassen ~ Top + Bottom, data = data.svm, kernel = "sigmoid") plot(fit2, data.svm, symbolPalette = mypalette[c(1,2)], svSymbol = "o", grid = 200, col = brewer.pal(3, "Pastel1")[c(3,1)]) ################################################### ### chunk number 23: ROCR ################################################### library(MASS) fit.lda <- lda(klassen ~ Bottom + Top, data = banknote.ausw) pred.lda <- predict(fit.lda, banknote.ausw) library(ROCR) pred <- prediction(pred.lda$x, klassen) ################################################### ### chunk number 24: ROCRplot1 ################################################### perf1 <- performance(pred, "tpr", "fpr") plot(perf1, main = "Klassische ROC-Kurve") abline(a = 0, b = 1) ################################################### ### chunk number 25: ROCRplot2 ################################################### perf2 <- performance(pred, "prec", "rec") plot(perf2, main = "Precision/recall Graph") ################################################### ### chunk number 26: AUC ################################################### performance(pred, "auc") ################################################### ### chunk number 27: ipred ################################################### library(ipred) mypredict.lda <- function(object, newdata) predict(object, newdata = newdata)$class errorest(klassen ~ Bottom + Top, data = banknote.ausw, model = lda, estimator = "cv", predict= mypredict.lda)