# STA 414, Spring 2014, Assignment #3 R script. source("pca.r") source("ass3-mlp.r") # READ THE DATA. trny <- scan("a3trny") trnx <- as.matrix(read.table("a3trnx",head=FALSE)) tsty <- scan("a3tsty") tstx <- as.matrix(read.table("a3tstx",head=FALSE)) # SPLIT TRAINING SET INTO ESTIMATION AND VALIDATION SETS. nest <- 1000 esty <- trny[1:nest] estx <- trnx[1:nest,] valy <- trny[-(1:nest)] valx <- trnx[-(1:nest),] # FIND PRINCIPAL COMPONENTS AND PROJECTIONS ON THEM. pc <- pca.vectors(estx,40) estpc <- pca.proj(pc,estx) valpc <- pca.proj(pc,valx) tstpc <- pca.proj(pc,tstx) # SET UP FOR PLOTTING TO A FILE. pdf("ass3-plot.pdf",height=7,width=5) par(mfrow=c(3,3),mar=c(4,3,4,0.2),mgp=c(1.3,0.5,0)) # SHOW SOME PLOTS OF PC'S (JUST TO SEE WHAT THEY'RE LIKE). plot(estpc[,1],estpc[,2],pch=20,col=c("blue","red")[esty+1]) title("Pairwise scatterplots of PCs",cex.main=0.8) plot(estpc[,3],estpc[,4],pch=20,col=c("blue","red")[esty+1]) plot(estpc[,39],estpc[,40],pch=20,col=c("blue","red")[esty+1]) # FUNCTION TO DISPLAY TEST ERROR RATE AND -LOG PROBABILITY. test_performance <- function (p) { # Set the condition to FALSE when manually tuning, so nothing will be # displayed that would let you "cheat". if (TRUE) { cat(" Test error rate =",round(mean((p>0.5)!=tsty),3), " -log pr =",round(mean(-log(p^tsty*(1-p)^(1-tsty))),3)) cat("\n") } } # FIT LOGISTIC REGRESSION MODELS WITH 10, 20, AND 40 PC'S. for (k in c(10,20,40)) { m <- glm(esty~estpc[,1:k],family="binomial") b <- coef(m) p <- 1 / (1 + exp( - (b[1] + valpc[,1:k]%*%b[-1]))) cat("\nLOGISTIC REGRESSION,", k, "PCS: validation error rate =", round(mean((p>0.5)!=valy),3), " -log pr =", round(mean(-log(p^valy*(1-p)^(1-valy))),3)) cat("\n") test_performance (1 / (1 + exp( - (b[1] + tstpc[,1:k]%*%b[-1])))) # Uncomment this to see the regression coefficients and other info. # print(summary(m)) } # FUNCTION TO COMPUTE AVERAGE LOG PROBABILITY FROM MLP ON VALIDATION CASES. valE <- function (fit, valx, valy) { E <- rep (NA, nrow(fit$W)) p <- ncol(valx) m <- (ncol(fit$W)-1) / (p+2) skel <- mlp_skeleton(p,m) for (i in 1:length(E)) { f <- mlp_forward (valx, relist(fit$W[i,],skel)) $ o E[i] <- mean (log (1+exp(-f*(2*valy-1)))) } E } # FIT NEURAL NETWORK MODELS WITHOUT PENALTY. m <- 10 # Number of hidden units for (k in c(10,20,40)) { set.seed(1) fit <- mlp_train(esty, estpc[,1:k], m, 0.00025, 30000) plot(fit$E,type="l",col="blue",ylim=c(0,0.5)) if (k==10) title("MLP, no penalty (k = 10, 20, 40)",cex.main=0.8) ve <- valE(fit,valpc[,1:k],valy) lines(ve,col="red") ve_min <- which.min(ve) wl <- relist (fit$W[ve_min,], mlp_skeleton(k,m)) p <- 1 / (1 + exp(-mlp_forward(valpc[,1:k],wl)$o)) cat("\nMLP,",k,"PCS:", "at",ve_min,"validation error rate =", round(mean((p>0.5)!=valy),3), " -log pr =", round(mean(-log(p^valy*(1-p)^(1-valy))),3)) cat("\n") test_performance (1 / (1 + exp(-mlp_forward(tstpc[,1:k],wl)$o))) } # FIT NEURAL NETWORK MODELS WITH PENALTY. Penalty magnitudes have been # manually adjusted to get what look like reasonable results. m <- 10 # Number of hidden units for (k in c(10,20,40)) { set.seed(1) lambda1 <- if (k==10) 0.15 else if (k==20) 10.0 else if (k==40) 0.4 fit <- mlp_train(esty, estpc[,1:k], m, 0.00025, 30000, lambda1) plot(fit$E,type="l",col="blue",ylim=c(0,0.5)) if (k==10) title("MLP, with penalty (k = 10, 20, 40)",cex.main=0.8) ve <- valE(fit,valpc[,1:k],valy) lines(ve,col="red") ve_min <- which.min(ve) wl <- relist (fit$W[ve_min,], mlp_skeleton(k,m)) p <- 1 / (1 + exp(-mlp_forward(valpc[,1:k],wl)$o)) cat("\nMLP,",k,"PCS,",lambda1,"PENALTY:", "at",ve_min,"validation err rate =", round(mean((p>0.5)!=valy),3), " -log pr =", round(mean(-log(p^valy*(1-p)^(1-valy))),3)) cat("\n") test_performance (1 / (1 + exp(-mlp_forward(tstpc[,1:k],wl)$o))) } dev.off()