# Classification drive <- "D:" code.dir <- paste(drive, "STAT5703/Code", sep = "/") data.dir <- paste(drive, "STAT5703/Data", sep = "/") # source(paste(code.dir, "makesyn.r", sep = "/")) source(paste(code.dir, "3DRotations.r", sep = "/")) source(paste(code.dir, "Pt_to_Array.r", sep = "/")) source(paste(code.dir, "CreateGaussians.r", sep = "/")) source(paste(code.dir, "CreateGrid.r", sep = "/")) source(paste(code.dir, "confusion_ex.r", sep = "/")) source(paste(code.dir, "example_display.r", sep = "/")) #======================== x <- 10 ex <- expression(x^2) eval(ex) example.display <- function(train.data, train.class, test.data, test.class, numb.xy.pts, title, model.exp, predict.exp, extra={}){ # Plot the data (coloured numbers) plot.class(train.data, train.class, test.data, test.class, main=title) numb.classes <- length(unique(train.class)) model <- plot.class.boundary(train.data, train.class, test.data, test.class, model.exp, predict.exp, numb.xy.pts, extra=extra) cat("===== The model is =====\n") print(model) cat("\n===== Training =====\n") points <- train.data pred.train.class <- eval(predict.exp) print(confusion.expand(pred.train.class, train.class)) # Print the total number of incorrect classifications misclass.train.ind <- which(pred.train.class!=train.class) cat("\n===== Misclassified training cases =====\n") if (length(misclass.train.ind) > 0) { tmp <- rbind(train.data[misclass.train.ind,]) rownames(tmp) <- misclass.train.ind # Print the incorrect classes print(tmp) # Display the incorrect classes in blue points(train.data[misclass.train.ind, ], col="blue", pch=24) } else { cat("No misclassified cases\n") } if (length(test.data) > 0) { cat("\n===== Test =====\n") points <- test.data pred.test.class <- eval(predict.exp) ind.tmp <- ((points-1) != test.class)*(1:length(test.class)) print(confusion.expand(pred.test.class, test.class)) misclass.test.ind <- which(pred.test.class!=test.class) cat("\n===== Misclassified test cases =====\n") if (length(misclass.test.ind) > 0) { tmp <- rbind(test.data[misclass.test.ind,]) rownames(tmp) <- misclass.test.ind # Print the incorrect classes print(tmp) # Display the incorrect classes in blue points(tmp, col="blue", pch=25) } else { cat("No misclassified cases\n") } } cat("===================================\n\n") bringToTop(which = dev.cur(), stay = FALSE) model } # Cloud 1 s.1 <- 500 c.1 <- Make.Gaussian(s.1, 0, 0, varx = 5, vary = 3, pi/4) c.1 <- add.arr.vec(c.1, c(3, 15)) # Cloud 2 s.2 <- 400 c.2 <- Make.Gaussian(s.2, 0, 0, varx = 5, vary = 3, pi/4) c.1 <- add.arr.vec(c.1, c(10, -5)) temp <- list(c.1 = c.1, s.1 = s.1, c.2 = c.2, s.2 = s.2) # save(temp, file="D:/DATA/Data Mining Data/cloud2distinct.rdata") load("D:/DATA/Data Mining Data/cloud2distinct.rdata") data.1 <- temp$c.1 sz.1 <- temp$s.1 data.2 <- temp$c.2 sz.2 <- temp$s.2 # Combine and plot G.1.data <- rbind(data.1, data.2) dimnames(G.1.data) <- list(NULL, c("X", "Y")) plot(G.1.data, col = c(rep(2, sz.1), rep(3, sz.2)), asp = 1, pch = 20) #================================== # Cloud 3 s.3 <- 500 c.3 <- Make.Gaussian(s.3, 0, 0, varx = 5, vary = 3, pi/4) c.3 <- add.arr.vec(c.3, c(3, 10)) # Cloud 4 s.4 <- 400 c.4 <- Make.Gaussian(s.4, 0, 0, varx = 5, vary = 3, pi/4) c.4 <- add.arr.vec(c.4, c(10, 3)) temp <- list(c.3 = c.3, s.3 = s.3, c.4 = c.4, s.4 = s.4) # save(temp, file="D:/DATA/Data Mining Data/cloud2lap.rdata") load("D:/DATA/Data Mining Data/cloud2lap.rdata") data.3 <- temp$c.3 sz.3 <- temp$s.3 data.4 <- temp$c.4 sz.4 <- temp$s.4 # Combine and plot G.2.data <- rbind(data.3, data.4) dimnames(G.2.data) <- list(NULL, c("X", "Y")) plot(G.2.data, col = c(rep(2, sz.3), rep(3, sz.4)), asp = 1, pch = 20) # Cloud 5 s.5 <- 600 c.5 <- Make.Gaussian(s.5, 0, 0, varx = 6, vary = 2, pi/3) c.5 <- add.arr.vec(c.5, c(-10,15)) temp <- list(c.5 = c.5, s.5 = s.5) # save(temp, file="D:/DATA/Data Mining Data/cloud5.rdata") load("D:/DATA/Data Mining Data/cloud5.rdata") sz.5 <- temp$s.5 data.5 <- temp$c.5 G.3.data <- rbind(data.3, data.4, data.5) dimnames(G.3.data) <- list(NULL, c("X", "Y")) # plot(G.3.data, col = c(rep(3, sz.3), rep(2, sz.4), rep(4, sz.5)), asp = 1, pch = 20) # 3 clouds C.1 <- Make.Gaussian(300, 8, 6, 4, 1, -pi/3) C.2 <- Make.Gaussian(500, -4, 6, 4, 1.2, pi/3) C.3 <- Make.Gaussian(300, 1, -15, 4, 1, 0) K1 <- rep(1, 300) K2 <- rep(2, 500) K3 <- rep(3, 300) tmp <- list(C.1 = C.1, C.2 = C.2, C.3 = C.3, K1 = K1, K2 = K2, K3 = K3) # save(tmp, file="D:/DATA/Data Mining Data/cloud3.rdata") load("D:/DATA/Data Mining Data/cloud3.rdata") data.6 <- tmp$C.1 data.7 <- tmp$C.2 data.8 <- tmp$C.3 sz.6 <- length(tmp$K1) sz.7 <- length(tmp$K2) sz.8 <- length(tmp$K3) data.123 <- data.frame(rbind(data.6, data.7, data.8)) plot(data.123, asp=1, col = c(rep(2, sz.6), rep(3, sz.7), rep(4, sz.8)), pch = 20) # Synthetic data.sz <- 500 res <- create.synthetic(data.sz) # Do the first time # save(res, file=paste(data.dir, "syn.Rdata", sep="/")) load(paste(data.dir, "syn.Rdata", sep="/")) # graphics.off() for (i in 1:8) { plot.class(res$data,res$tp[,i], {}, {}) readline("Press Enter...") } #========= Nearest Neighbors ================== # Set the number of neighbours wanted f.create.data <- function(numb.classes, n.in.class, centre.class) { total.n.in.classes <- sum(n.in.class) data <- matrix(0, total.n.in.classes, 3) start <- 1 for (i in (1:numb.classes)) { end <- start + n.in.class[i] -1 pos <- start:end data[pos,1] <- rnorm(n.in.class[i], 0, 1) + centre.class[[i]][1] data[pos,2] <- rnorm(n.in.class[i], 0, 1) + centre.class[[i]][2] data[pos,3] <- rep(i, n.in.class[i]) start <- end + 1 } data } draw.circle <- function(r, x0, y0, c) { t <- seq(0,2*pi,by=.01) x <- x0 + r*cos(t) y <- y0 + r*sin(t) lines (x,y, col=c) } data <- f.create.data(3, c(50, 60, 70), list(c(3,3), c(5,5), c(5,3))) grid <- f.create.grid(data, c(10,10)) # myknn <- function(k, data, grid) { a.min <- min(min(grid[[1]]), min(grid[[2]])) a.max <- max(max(grid[[1]]), max(grid[[2]])) plot(data[,1], data[,2], xlim=c(a.min, a.max), ylim=c(a.min, a.max), col = data[,3]) # dist <- matrix(0, 2, k); numb.classes <- length(unique(data[,3])) curr.class <- matrix(0, 1, numb.classes) p.class <- matrix(0, 1, dim(grid[[3]])[[1]]) # # Draw a circle in white - doesn't show. # Set the x, y, and radius i.e. max.old # draw.circle(0, 0, 0, 'white') x.old <- 0; y.old <- 0; max.old <- 0 # for (j in (1:(dim(grid[[3]])[[1]]))){ ## run through all the test points grid for (i in 1: dim(data)[[1]]) { currdist <- (data[i,1]-grid[[3]][j,1])*(data[i,1]-grid[[3]][j,1])+(data[i,2]-grid[[3]][j,2])*(data[i,2]-grid[[3]][j,2]); if (i <= k) { # Set the first k dist[1,i] <- currdist; dist[2,i] <- i; } # if (i <= k) else {## go through the k list and replace the biggest one if possible max.d <- max(dist[1,]) max.index <- which.max(dist[1,]) if (currdist < max.d) { dist[1,max.index] <- currdist; dist[2,max.index] <- i; } } } for (cl in (1:numb.classes)) { curr.class[cl] <- sum(data[,3][dist[2,]]== cl) } # take a vote v.max <- curr.class==max(curr.class) v <- sum(v.max) # v = 1 no tie v.index <- which(v.max) p.class[j] <- v.index[floor(runif(1)/v)+1] draw.circle(max.old^(.5), x.old, y.old, 'white') x.old <- grid[[3]][j,1] y.old <- grid[[3]][j,2] max.old <- max.d draw.circle(max.d^(.5), grid[[3]][j,1], grid[[3]][j,2], 'black') points(data[,1], data[,2], col = data[,3]) points(grid[[3]][,1], grid[[3]][,2], col = p.class, pch = 20) if (as.numeric(rownames(grid$zp)[j])%%14 == 5) { readline("Press return") } } contour(grid[[1]],grid[[2]],matrix(p.class,nrow=length(grid[[1]])), levels = 1:(numb.classes-1)+.5) points(data[,1], data[,2], col = data[,3]) points(grid[[3]][,1], grid[[3]][,2], col = p.class, pch=20) p.class } graphics.off() my.class <- myknn(5, data, grid) #========================== # Nearest neighbour #========================== (nn.3 <- knn(data[ ,1:2], cbind(10:80/10, rep(3.5,71)), data[ ,3], k = 3)) unclass(nn.3) unclass(nn.3)[1]+2 ################################################# # k-nearest neighbour classification for test set # from training set. For each row of the test set # the k nearest (in Euclidean distance) training # set vectors are found, and the classification is # decided by majority vote, with ties broken at random. # If there are ties for the kth nearest vector, # all candidates are included in the vote. ################################################# xa <- c(-2,-1,1,2) ya <- c(-2,-1,1,2) class <- c(0,0,1,1,0,0,1,1,1,1,0,0,1,1,0,0) # model.exp <- expression({}) predict.exp <- expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1) graphics.off() oldpar <- par(mfrow=c(2,3)) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "1-NN", model.exp, predict.exp, 1) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "2-NN", model.exp, predict.exp, 2) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "3-NN", model.exp, predict.exp, 3) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "4-NN", model.exp, predict.exp, 4) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "5-NN", model.exp, predict.exp, 5) example.display(expand.grid(xa, ya), class, {}, {}, c(100, 100), "6-NN", model.exp, predict.exp, 6) par(oldpar) #============================ # 4 classes xb <- c(-2,-1,1,2) yb <- c(-2,-1,1,2) class.b <- c(0,0,2,2,0,0,2,2,1,1,3,3,1,1,3,3) # graphics.off() x11(height = 4, width = 8) oldpar <- par(mfrow=c(1,3)) example.display(expand.grid(xb, yb), class.b, {}, {}, c(100, 100), "2-NN", model.exp, predict.exp, 2) example.display(expand.grid(xb, yb), class.b, {}, {}, c(100, 100), "3-NN", model.exp, predict.exp, 3) example.display(expand.grid(xb, yb), class.b, {}, {}, c(100, 100), "4-NN", model.exp, predict.exp, 4) par(oldpar) #============================================ # A function to produce those indices NOT in a set. # Use this to get the test sample. #============================================ "%w/o%" <- function(x,y) x[!x %in% y] #============================================ # Set the indices for the training/test sets #============================================ get.train <- function (data.sz, train.sz) { # Take subsets of data for training/test samples # Return the indices train.ind <- sample(data.sz, train.sz) test.ind <- (1:data.sz) %w/o% train.ind list(train=train.ind, test=test.ind) } tt.ind <- get.train(data.sz, 400) # save(tt.ind, file=paste(data.dir, "synTrainTest.Rdata", sep="/")) load(paste(data.dir, "synTrainTest.Rdata", sep="/")) syn.train.class <- res$tp[tt.ind[[1]],] syn.train.data <- res$data[tt.ind[[1]],] # syn.test.class <- res$tp[tt.ind[[2]],] syn.test.data <- res$data[tt.ind[[2]],] f.menu <- function(){ while (TRUE) { cat("Enter the number of the data set","\n") cat ("0 to quit\n") resN <- menu(c("1:2", "1/3:2/4", "1:2/3:4", "1:2 rotated", "1/3:2/4 rotated", "1:2/3:4 rotated", "Hole", "x^2")) if (resN == 0) break size <- readline("# of neighbours - 1 to 10: ") example.display(syn.train.data, syn.train.class[ ,resN] , {}, {}, c(100, 100), paste(as.integer(size),"- NN"), expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), as.integer(size)) } } f.menu() #====================== End KNN =========== # Logistic regression p.hat <- function(object, X) { X <- as.matrix(X) term <- cbind(rep(1,dim(X)[1]), X)%*%as.matrix(object$coefficients) 1/(1+exp(-term)) } (log.1 <- glm(c(rep(0, sz.1), rep(1, sz.2)) ~ G.1.data[,1] + G.1.data[,2], quasibinomial(link = "logit"))) summary(log.1) plot(G.1.data, col = c(rep(2, sz.1), rep(3, sz.2)), asp = 1, pch = 20) abline(-log.1[[1]][1]/log.1[[1]][3], -log.1[[1]][2]/log.1[[1]][3]) # model.exp <- expression(glm(factor(train.class) ~ train.data[,1] + train.data[,2], quasibinomial(link = "logit"))) predict.exp <- expression((p.hat(model, points) > 0.5)*1) lr.1 <- example.display(G.1.data[,1:2], c(rep(0,sz.1), rep(1,sz.2)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) #================================== lr.1$coefficients #================================== lr.1$df.null #================================== lr.1$df.residual #================================== lr.1$null.deviance #================================== lr.1$deviance example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) #============================ # First 3 clouds #=========================== (LR3.1 <- example.display(G.3.data[,1:2], c(rep(1,sz.3), rep(1,sz.4), rep(0,sz.5)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp)) #=========================== (LR3.2 <- example.display(G.3.data[,1:2], c(rep(1,sz.3), rep(0,sz.4), rep(1,sz.5)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp)) #=========================== plot(G.3.data, col = c(rep(3, sz.3), rep(2, sz.4), rep(4, sz.5)), asp = 1, pch = 20) abline(-LR3.1$coeff[1]/LR3.1$coeff[3], -LR3.1$coeff[2]/LR3.1$coeff[3]) abline(-LR3.2$coeff[1]/LR3.2$coeff[3], -LR3.2$coeff[2]/LR3.2$coeff[3]) # nn example.display(G.3.data[,1:2], c(rep(0, sz.3), rep(1, sz.4), rep(2, sz.5)), {}, {}, c(100, 100), "1 - NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = 1)) - 1)) # Second 3 cloud predict.exp <- expression((p.hat(model, points) > 0.5)*1) LR3.3 <- example.display(data.123, c(rep(0, sz.6), rep(1, sz.7), rep(1, sz.8)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) predict.exp <- expression((p.hat(model, points) > 0.5)*2) LR3.4 <- example.display(data.123, c(rep(0, sz.6), rep(0, sz.7), rep(2, sz.8)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) predict.exp <- expression((p.hat(model, points) > 0.5)*1 + 1) LR3.5 <- example.display(data.123, c(rep(2, sz.6), rep(1, sz.7), rep(2, sz.8)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) plot(data.123, col = c(rep(3, sz.6), rep(2, sz.7), rep(4, sz.8)), asp = 1, pch = 20) abline(-LR3.3$coeff[1]/LR3.3$coeff[3], -LR3.3$coeff[2]/LR3.3$coeff[3]) abline(-LR3.4$coeff[1]/LR3.4$coeff[3], -LR3.4$coeff[2]/LR3.4$coeff[3]) abline(-LR3.5$coeff[1]/LR3.5$coeff[3], -LR3.5$coeff[2]/LR3.5$coeff[3]) ############# NN ############### example.display(data.123, c(rep(0, sz.6), rep(1, sz.7), rep(2, sz.8)), {}, {}, c(100, 100), "11 - NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = 11)) - 1)) ######### pairwise ############## predict.exp <- expression((p.hat(model, points) > 0.5)*1) temp <- data.frame(rbind(data.6, data.7)) mod.67 <- example.display(temp, c(rep(0, sz.6), rep(1, sz.7)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) predict.exp <- expression((p.hat(model, points) > 0.5)*1 + 1) temp <- data.frame(rbind(data.7, data.8)) mod.78 <- example.display(temp, c(rep(1, sz.7), rep(2, sz.8)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) predict.exp <- expression((p.hat(model, points) > 0.5)*2) temp <- data.frame(rbind(data.8, data.6)) mod.86 <- example.display(temp, c(rep(2, sz.8), rep(0, sz.6)), {}, {}, c(100,100), "Logistic Regression", model.exp, predict.exp) res <- f.create.grid(data.123, c(100,100)) xp <- res$xp yp <- res$yp points <- res$zp dimnames(points)[[2]]<- dimnames(data.123)[[2]] comb.class <- matrix(0, dim(points)[1], 3) Z <- (p.hat(mod.67, points) > 0.5)*1 comb.class[which(Z == 0), 1] <- comb.class[which(Z == 0), 1] + 1 comb.class[which(Z == 1), 2] <- comb.class[which(Z == 1), 2] + 1 Z <- (p.hat(mod.78, points) > 0.5)*1 + 1 comb.class[which(Z == 1), 2] <- comb.class[which(Z == 1), 2] + 1 comb.class[which(Z == 2), 3] <- comb.class[which(Z == 2), 3] + 1 Z <- (p.hat(mod.86, points) > 0.5)*2 comb.class[which(Z == 2), 3] <- comb.class[which(Z == 2), 3] + 1 comb.class[which(Z == 0), 1] <- comb.class[which(Z == 0), 1] + 1 plot(data.123, asp=1, col = c(rep(2, sz.6), rep(3, sz.7), rep(4, sz.8)), pch = 20) points(points, col = max.col(comb.class) + 1, pch = ".", cex = 2) abline(-mod.67$coeff[1]/mod.67$coeff[3], -mod.67$coeff[2]/mod.67$coeff[3]) abline(-mod.78$coeff[1]/mod.78$coeff[3], -mod.78$coeff[2]/mod.78$coeff[3]) abline(-mod.86$coeff[1]/mod.86$coeff[3], -mod.86$coeff[2]/mod.86$coeff[3]) #======= NN vs Logistic==================== # Compare graphics.off() oldpar <- par(mfrow=c(2,2), mar = c(2, 2, 2, 1)) example.display(G.1.data[,1:2], c(rep(0,sz.1), rep(1,sz.2)), {}, {}, c(100,100), "Logistic Regression", expression(glm(factor(train.class)~train.data[,1]+train.data[,2], quasibinomial(link = "logit"))), expression((p.hat(model, points) > 0.5)*1)) example.display(G.1.data[,1:2], c(rep(0,sz.1), rep(1,sz.2)), {}, {}, c(100, 100), "3-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 3) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100,100), "Logistic Regression", expression(glm(factor(train.class)~train.data[,1]+train.data[,2], quasibinomial(link = "logit"))), expression((p.hat(model, points) > 0.5)*1)) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100, 100), "3-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 3) par(oldpar) graphics.off() oldpar <- par(mfrow=c(2,2), mar = c(2, 2, 2, 1)) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100, 100), "5-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 5) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100, 100), "7-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 7) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100, 100), "9-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 9) example.display(G.2.data[,1:2], c(rep(0,sz.3), rep(1,sz.4)), {}, {}, c(100, 100), "11-NN", expression({}), expression(unclass(knn(train.data, data.frame(points), train.class, k = extra)) - 1), 11) par(oldpar) f.menu <- function(){ while (TRUE) { cat("Enter the number of the data set","\n") cat ("0 to quit\n") resN <- menu(c("1:2", "1/3:2/4", "1:2/3:4", "1:2 rotated", "1/3:2/4 rotated", "1:2/3:4 rotated", "Hole", "x^2")) if (resN == 0) break example.display(syn.train.data, syn.train.class[,resN], syn.test.data, syn.test.class[,resN], c(100, 100), "Logistic Regression", model.exp, predict.exp) } } f.menu() #====================== LDA/QDA =========== source(paste(code.dir, "creategaussians.r", sep="/")) source(paste(code.dir, "3Drotations.r", sep="/")) library(MASS) cols <- c(1,2,3) C1 <- Make.Gaussian(300, 1, 6, 4, 1, pi/3) C2 <- Make.Gaussian(500, 1, -6, 4, 1.2, pi/3) C3 <- Make.Gaussian(300, 1, 6, 2, 1, pi/3) C4 <- Make.Gaussian(500, 1, -6, 8, 1, pi/3) cov.diff <- list(X1 = C1, X2 = C2, X3 = C3, X4 = C4) # save(cov.diff, file="D:/DATA/Data Mining Data/ldaqdacovdiff.rdata") load("C:/Documents and Settings/Kathryn/My Documents/Courses/STAT5703/Data/ldaqdacovdiff.rdata") X.1 <- cov.diff$X1 n.1 <- length(X.1[ ,1]) colnames(X.1) <- c("x1", "x2") X.2 <- cov.diff$X2 n.2 <- length(X.2[ ,1]) colnames(X.2) <- c("x1", "x2") C.1 <- rep(1,300) C.2 <- rep(2,500) X.12 <- data.frame(rbind(X.1, X.2)) (sig.1 <- var(X.1)) (sig.2 <- var(X.2)) (sig.12 <- ( (n.1 - 1)*sig.1 + (n.2 - 1)*sig.2 ) / (n.1 + n.2 - 2)) C.12 <- c(C.1, C.2) plot(X.12, asp=1, col = C.12 + 1, pch = 20) X.3 <- cov.diff$X3 n.3 <- length(X.3[,1]) colnames(X.3) <- c("x1", "x2") X.4 <- cov.diff$X4 n.4 <- length(X.4[,1]) colnames(X.4) <- c("x1", "x2") C.3 <- rep(1, 300) C.4 <- rep(2, 500) X.34 <- data.frame(rbind(X.3, X.4)) (sig.3 <- var(X.3)) (sig.4 <- var(X.4)) (sig.34 <- ( (n.3 - 1)*sig.3 + (n.4 - 1)*sig.4 ) / (n.3 + n.4 - 2)) C.34 <- c(C.3, C.4) plot(X.34, asp=1, col = C.34 + 1, pch = 20) # Group means (mu.1 <- matrix(apply(X.1, 2, mean))) (mu.2 <- matrix(apply(X.2, 2, mean))) (mu.3 <- matrix(apply(X.3, 2, mean))) (mu.4 <- matrix(apply(X.4, 2, mean))) # Covariance from scratch ((t(X.1) - drop(mu.1))%*%t(t(X.1) - drop(mu.1)))/(n.1 - 1) (mu.diff.12 <- mu.1 - mu.2) (mu.diff.34 <- mu.3 - mu.4) (tmp <- svd(sig.12)) tmp$v %*% diag(1/tmp$d) %*% t(tmp$u) (sig.inv.12 <- ginv(sig.12)) (sig.inv.34 <- ginv(sig.34)) (l.12 <- sig.inv.12 %*% (mu.diff.12)) (l.34 <- sig.inv.34 %*% (mu.diff.34)) Y.12 <- t(l.12)%*% t(as.matrix(X.12)) Y.34 <- t(l.34)%*% t(as.matrix(X.34)) plot(Y.12, 1:800, col = C.12 + 1, pch = 20) rug(Y.12, col = 4) plot(Y.34, 1:800, col = C.12 + 1, pch = 20) rug(Y.34, col = 4) mid.12 <- (mu.1 + mu.2)/2 (m.12 <- 0.5*drop( t(l.12)%*%as.matrix(mid.12)) ) mid.34 <- (mu.3 + mu.4)/2 (m.34 <- 0.5*drop( t(l.34)%*%as.matrix(mid.34)) ) plot(Y.12, 1:800, col = C.12 + 1, pch = 20) lines(c(m.12, m.12), c(0, 800), lwd = 2, col = 4) rug(Y.12, col = 4) (X.12.lda <- lda(X.12, C.12)) mu.1 mu.2 n.1/(n.1 + n.2) n.2/(n.1 + n.2) (sl.12 <- l.12[2,1]/l.12[1,1]) (sl.12.lda <- X.12.lda$scaling[2,1]/X.12.lda$scaling[1,1]) X.G <- seq(-15, 15, length = 100) Y.G <- X.G Z.G <- expand.grid(X.G, Y.G) colnames(Z.G) <- colnames(X.12) Z.G.cl <- as.numeric(predict(X.12.lda, Z.G)$class) plot(Z.G, col = cols[Z.G.cl] + 1, xlim = c(-15, 15), ylim = c(-15, 15), xlab = expression(x[1]), ylab = expression(x[2]), pch = ".", main = "LDA", asp = 1) points(X.12, col = C.12 + 1, pch = 20) lines(-15:15, sl.12*(-15:15) - 5, lwd = 2) lines(-9:9, -1/sl.12*(-9:9 - mid.12[1,1]) + mid.12[2,1], , lwd = 2, col = "blue") (X.34.lda <- lda(X.34, C.34)) (sl.34 <- l.34[2,1]/l.34[1,1]) (sl.34.lda <- X.34.lda$scaling[2,1]/X.34.lda$scaling[1,1]) Z.G.cl <- as.numeric(predict(X.34.lda, Z.G)$class) plot(Z.G, col = cols[Z.G.cl] + 1, xlim = c(-15, 15), ylim = c(-15, 15), xlab = expression(x[1]), ylab = expression(x[2]), pch = ".", main = "LDA", asp = 1) points(X.34, col = C.34 + 1, pch = 20) lines(-15:15, sl.34*(-15:15) - 5, lwd = 2) lines(-9:9, -1/sl.34*(-9:9 - mid.34[1,1]) + mid.34[2,1], , lwd = 2, col = "blue") # QDA (sig.inv.1 <- ginv(sig.1)) (sig.inv.2 <- ginv(sig.2)) (sig.inv.3 <- ginv(sig.3)) (sig.inv.4 <- ginv(sig.4)) (t(mu.1) %*% sig.inv.1 - t(mu.2) %*% sig.inv.2) 0.5*(sig.inv.1 - sig.inv.2) 0.5*(t(mu.1) %*% sig.inv.1 %*% mu.1 - t(mu.2) %*% sig.inv.2 %*% mu.2) 0.5*log(det(sig.1)/det(sig.2)) (X.12.qda <- qda(X.12, C.12)) X.G <- seq(-15, 15, length = 100) Y.G <- X.G Z.G <- expand.grid(X.G, Y.G) colnames(Z.G) <- colnames(X.12) Z.G.cl <- as.numeric(predict(X.12.qda, Z.G)$class) plot(Z.G, col = cols[Z.G.cl] + 1, xlim = c(-15, 15), ylim = c(-15, 15), xlab = expression(x[1]), ylab = expression(x[2]), pch = ".", main = "QDA", asp = 1) points(X.12, col = C.12 + 1, pch = 20) # Plot the curve x.1 <- seq(-15, 15, 0.1) x.2 <- 2.735*x.1+145.6-5.4148*10^(-4)*sqrt(3.4064*10^(6)*x.1^2+9.6527*10^(8)*x.1+7.1404*10^(10)) lines(x.1, x.2, col = "blue", lwd = 2) # QDA on 34 (t(mu.3) %*% sig.inv.3 - t(mu.4) %*% sig.inv.4) 0.5*(sig.inv.3 - sig.inv.4) 0.5*(t(mu.3) %*% sig.inv.3 %*% mu.3 - t(mu.4) %*% sig.inv.4 %*% mu.4) 0.5*log(det(sig.3)/det(sig.4)) (X.34.qda <- qda(X.34, C.34)) Z.G.cl <- as.numeric(predict(X.34.qda, Z.G)$class) plot(Z.G, col = cols[Z.G.cl] + 1, xlim = c(-15, 15), ylim = c(-15, 15), xlab = expression(x[1]), ylab = expression(x[2]), pch = ".", main = "QDA", asp = 1) points(X.34, col = C.34 + 1, pch = 20) x <- seq(-14, 6, 0.1) y.2 <- (-41102*x+2954900-2*sqrt( -894656076*x^2-264075242400*x+2192040941305))/79823 lines(x,y.2, col = "blue", lwd = 2) # Angle X.5 <- Make.Gaussian(300, -6, 36, 2, 1, pi/4) X.6 <- Make.Gaussian(500, 1, -5, 8, 1, -pi/4) X.56 <- data.frame(rbind(X.5, X.6)) predict.exp <- expression(unclass(predict(model, points)$class)-1) qda.56 <- example.display(X.56, c(rep(0, 300), rep(1, 500)), {}, {}, c(100, 100), "QDA", expression(qda(train.data, train.class)), predict.exp) predict.exp <- expression(unclass(predict(model, points)$class)-1) lda.1 <- example.display(data.123, c(rep(0, sz.6), rep(1, sz.7), rep(2, sz.8)), {}, {}, c(100, 100), "LDA", expression(lda(train.data, train.class)), predict.exp) predict(lda.1, cbind(-15:15, rep(-10, 31))) qda.1 <- example.display(data.123, c(rep(0, sz.6), rep(1, sz.7), rep(2, sz.8)), {}, {}, c(100, 100), "QDA", expression(qda(train.data, train.class)), predict.exp) predict(qda.1, cbind(-15:15, rep(-10, 31))) (data.123.qda <- qda(data.123, c(rep(2, sz.6), rep(3, sz.7), rep(4, sz.8)))) plot(data.123, asp=1, col = c(rep(2, sz.6), rep(3, sz.7), rep(4, sz.8)), pch = 20) points(points, col = max.col(comb.class) + 1, pch = ".", cex = 2) ################ predict.exp <- expression(unclass(predict(model, points)$class)-1) f.menu <- function(){ while (TRUE) { cat("Enter the type of classification - 0 to stop","\n") res <- menu(c("LDA", "QDA" )) if (res == 0) break cat("Enter the number of the data set","\n") resN <- menu(c("1:2", "1/3:2/4", "1:2/3:4", "1:2 rotated", "1/3:2/4 rotated", "1:2/3:4 rotated", "Hole", "1 rotated")) switch(res, example.display(syn.train.data, syn.train.class[,resN], syn.test.data, syn.test.class[,resN], c(100, 100), "LDA", expression(lda(train.data, train.class)), predict.exp), example.display(syn.train.data, syn.train.class[,resN], syn.test.data, syn.test.class[,resN], c(100, 100), "QDA", expression(qda(train.data, train.class)), predict.exp)) } } f.menu() #============= Fleas =============== source(paste(code.dir, "ReadFleas.r", sep="/")) flea.lda <- example.display(data.frame(d.flea[,6],d.flea[,1]), species-1, {}, {}, c(100, 100), "LDA", expression(lda(train.data, train.class)), predict.exp) ############################ flea.qda <- example.display(data.frame(d.flea[,6],d.flea[,1]), species-1, {}, {}, c(100, 100), "QDA", expression(qda(train.data, train.class)), predict.exp) # (d.diff <- matrix(c(c(1,1,1,1,1,2,2,2,2,2), c(1,1,1,1,2,2,2,2,2,2), c(1,1,1,2,2,2,2,2,2,2), c(1,1,1,2,2,3,2,2,2,2), c(1,1,2,2,2,3,3,2,2,2), c(1,1,2,2,2,3,3,3,2,2), c(1,2,2,2,2,3,3,3,3,2), c(1,2,2,2,3,3,3,3,3,2), c(2,2,2,2,3,3,3,3,3,2), c(2,2,2,2,3,3,3,3,3,3)), 10, 10, byrow=T)) plot(expand.grid(1:10, 1:10), col=d.diff[,10:1]+1, pch=16) (change.a <- apply(d.diff, 2, diff)) points(expand.grid(1:9 +0.5, 1:10), col=abs(change.a[,10:1]), pch=3) (change.b <- t(apply(d.diff, 1, diff))) points(expand.grid(1:10, 1:9+0.5), col=abs((change.b)[,9:1]), pch=3) disp.3d <- function(model, data, model.type, steps=30) { # Get a bounding interval for data[,3] # This will be taken as our third dimension z.min <- floor(min(data[,3])) z.max <- ceiling(max(data[,3])) # Break the interval into a number of subintervals z.int <- seq(z.min, z.max, by=((z.max-z.min)/steps)) # For the other two variables we create a grid res <- f.create.grid(data.frame(data[,1],data[,2]), c(steps,steps)) xp <- res$xp yp <- res$yp gp <- res$zp # The grid points # In order to display the separators we look for # changes in predicted class x.change <- {} y.change <- {} z.change <- {} for (i in 1:steps){ # We will look at a grids on a sequence of planes # in the d.flea[,4] direction. data.3d <- cbind(gp, z.int[i]) # The predict method likes the names to be the same colnames(data.3d) <- colnames(data) # We use the model to predict the class of every point # on the grid on the current data[,3] plane if (model.type == "lda") { Z <- unclass(predict(model, data.3d)$class) } else if (model.type == "qda") { Z <- unclass(predict(model, data.3d)$class) } else if (model.type == "net") { Z <- max.col(predict(model, data.3d)) } # Put the result into a matrix Z.m <- matrix(Z, nrow=steps+1) # Find the positions in a row where the class changes Z.m.diff.r <- t(apply(Z.m, 1, diff)) cc <- which(Z.m.diff.r != 0, arr.ind=T) # Find the positions in a column where the class changes Z.m.diff.c <- apply(Z.m, 2, diff) rr <- which(Z.m.diff.c != 0, arr.ind=T) # Convert the matrix positions to grid positions x.change <- c(x.change, xp[rr[,1]],xp[cc[,1]]) y.change <- c(y.change, yp[rr[,2]],yp[cc[,2]]) z.change <- c(z.change, rep(z.int[i],length(rr[,1])+length(cc[,1]))) } # After we have found the change points on all the grids data.frame(x.change, y.change, z.change) } data.for.3d <- data.frame(d.flea[,6],d.flea[,1],d.flea[,4]) colnames(data.for.3d) <- c("aede3","tars1","aede1") # Do an LDA on the three dimensional data (flea.lda <- lda(data.for.3d, species-1)) # Get the change points data.3d <- disp.3d(flea.lda, data.for.3d, "lda", steps=60) library(rgl) plot3d(d.flea[,c(6,1,4)], type="s", col = species+1, size=0.3) points3d(data.3d, size=2) (flea.qda <- qda(data.for.3d, species-1)) data.3d <- disp.3d(flea.qda, data.for.3d, "qda", steps=60) plot3d(d.flea[,c(6,1,4)], type="s", col = species+1, size=0.3) points3d(data.3d, size=3)