# STA 414/2104, Spring 2006, Assignment #2. # # R. M. Neal # # Uses the PCA and LDA functions from pca.r and lda.r. Note: The LDA function # handles any number of classes, and for two classes, uses 1 and 2 to identify # the classes, rather than 0 and 1, as in the data for this assignment. source("pca.r") source("lda.r") # Read the data. x.train <- as.matrix(read.table("a2-x-train",head=F)) x.test <- as.matrix(read.table("a2-x-test",head=F)) y.train <- scan("a2-y-train") y.test <- scan("a2-y-test") # Find first ten principal components. pc <- pca.vectors (x.train, 10) pc.train <- pca.proj (pc, x.train) pc.test <- pca.proj (pc, x.test) # Find block averages. ave.train <- matrix(0,nrow(x.train),ncol(x.train)/20) ave.test <- matrix(0,nrow(x.test),ncol(x.train)/20) for (j in 1:ncol(ave.train)) { ave.train[,j] <- apply(x.train[,(1+20*(j-1)):(20*j)],1,mean) ave.test[,j] <- apply(x.test[,(1+20*(j-1)):(20*j)],1,mean) } # Try logistic regression models. lr.ave <- glm (y.train ~ ave.train, family=binomial) lr.b.ave <- coef(lr.ave) p.lr.ave <- as.numeric (lr.b.ave[1] + ave.test%*%lr.b.ave[2:length(lr.b.ave)] > 0) lr.pc <- glm (y.train ~ pc.train[,1:10], family=binomial) lr.b.pc <- coef(lr.pc) p.lr.pc <- as.numeric (lr.b.pc[1] + pc.test[,1:10]%*%lr.b.pc[2:length(lr.b.pc)] > 0) # Try LDA. lda.ave <- lda.train (ave.train, y.train+1) p.lda.ave <- lda.predict (lda.ave, ave.test) - 1 lda.pc <- lda.train (pc.train[,1:10,drop=F], y.train+1) p.lda.pc <- lda.predict (lda.pc, pc.test[,1:10,drop=F]) - 1 # Print results. res <- c (lr.ave = mean(p.lr.ave!=y.test), lr.pc = mean(p.lr.pc!=y.test), lda.ave = mean(p.lda.ave!=y.test), lda.pc = mean(p.lda.pc!=y.test) ) print(res) # Put the next four plots on one page. par(mfrow=c(2,2)) # Compare the discriminant functions found by LDA and logistic regression # for block averages and principal components, in terms of the 10 variables # found using PCA or block averaging. plot (0:10, rep(0,11), ylim=range(c(lda.ave$a[2]-lda.ave$a[1],lda.ave$b[2,]-lda.ave$b[1,],lr.b.ave)), xlab="Block average discriminant functions for LDA (blue) and LR (green)", ylab="", type="n") lines(0:10, c(lda.ave$a[2]-lda.ave$a[1],lda.ave$b[2,]-lda.ave$b[1,]), col="blue") lines(0:10, lr.b.ave, col="green") abline(h=0) title("Reduced to 10 dimensions") plot (0:10, rep(0,11), ylim=range(c(lda.pc$a[2]-lda.pc$a[1],lda.pc$b[2,]-lda.pc$b[1,],lr.b.pc)), xlab="PC discriminant functions for LDA (blue) and LR (green)", ylab="", type="n") lines(0:10, c(lda.pc$a[2]-lda.pc$a[1],lda.pc$b[2,]-lda.pc$b[1,]), col="blue") lines(0:10, lr.b.pc, col="green") abline(h=0) title("Reduced to 10 dimensions") # Compare the discriminant functions found by LDA and logistic regression, # in terms of the 200 original inputs. oi.lda.ave <- rep (lda.ave$b[2,]-lda.ave$b[1,], each=20) oi.lr.ave <- rep (lr.b.ave[2:11], each=20) oi.lda.pc <- as.vector (pc$vectors %*% (lda.pc$b[2,]-lda.pc$b[1,])) oi.lr.pc <- as.vector (pc$vectors %*% lr.b.pc[2:11]) plot (0:200, rep(0,201), ylim=range(c(lda.ave$a[2]-lda.ave$a[1],oi.lda.ave,lr.b.ave[1],oi.lr.ave)), xlab="Block average discriminant functions for LDA (blue) and LR (green)", ylab="", type="n") abline(h=0) abline(v=seq(0.5,200.5,length=11),col="gray") lines(0:200, c(lda.ave$a[2]-lda.ave$a[1],oi.lda.ave),col="blue") lines(0:200, c(lr.b.ave[1],oi.lr.ave), col="green") title("Original 200 inputs") plot (0:200, rep(0,201), ylim=range(c(lda.pc$a[2]-lda.pc$a[1],oi.lda.pc,lr.b.pc[1],oi.lr.pc)), xlab="PC discriminant functions for LDA (blue) and LR (green)", ylab="", type="n") abline(h=0) abline(v=seq(0.5,200.5,length=11),col="gray") lines(0:200, c(lda.pc$a[2]-lda.pc$a[1],oi.lda.pc),col="blue") lines(0:200, c(lr.b.pc[1],oi.lr.pc),col="green") title("Original 200 inputs")