# STA 414, SPRING 2011, ASSIGNMENT 2 SOLUTION, R SCRIPT.
#
# Radford M. Neal, 2011

# Define the MLP functions, as modified for this assignent.

source("a2-mlp.txt")

# Read all the data, and put in suitable matrices.

est <- as.matrix(read.table("est.txt",head=F))
val <- as.matrix(read.table("val.txt",head=F))
tst <- as.matrix(read.table("tst.txt",head=F))

trn <- rbind(val,est)
trn.X <- trn[,-51]
trn.y <- trn[,51]

trn.mean <- apply(trn.X,2,mean)
trn.sd <- apply(trn.X,2,sd)
trn.X <- scale(trn.X,center=trn.mean,scale=trn.sd)

tst.y <- tst[,51]
tst.X <- tst[,-51]
tst.X <- scale(tst.X,center=trn.mean,scale=trn.sd)

# Print log probability of test cases, given predicted probabilities of a 1.

test.log.prob <- function (tst.y, p1)
{ 
  cat("Test log prob =", round(mean(ifelse(tst.y==1,log(p1),log(1-p1))),3),"\n")
}

# Try linear logistic regression.

lr <- glm (trn.y ~ trn.X, family="binomial")
lr.p1 <- 1 / (1 + exp(-(coef(lr)[1] + tst.X %*% coef(lr)[-1])))
cat("Logistic regression:  "); test.log.prob (tst.y, lr.p1)

# Plots go to PDF file, 8 per page.

pdf("a2-plots.pdf",pointsize=9,height=11,width=8)
par(mfrow=c(4,2))

title.string <- function (M,lambda1,lambda2,seed)
{ paste("M=",M,", lambda1=",lambda1,", lambda2=",lambda2,", seed=",seed,sep="")
}

# Train MLP networks with no penalty.

seeds <- 1:2
lambda1 <- lambda2 <- 0

for (M in c(2,4,8,16))
{ for (seed in seeds)
  { set.seed(seed)
    cat("\n")
    cat(title.string(M,lambda1,lambda2,seed),"\n")
    r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda1, lambda2, 1:nrow(val), 
                 3000, M, cv.plot=TRUE)
    title(title.string(M,lambda1,lambda2,seed))
    test.log.prob (tst.y, mlp.predict(tst.X,r))
  }
}

# Train MLP networks with 16 hidden units, with varying penalties (same for 
# both groups of weights).

M <- 16
seeds <- 1:2

for (lambda in c(1,3,9,27))
{ for (seed in seeds)
  { set.seed(seed)
    cat("\n")
    cat(title.string(M,lambda,lambda,seed),"\n")
    r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda, lambda, 1:nrow(val), 
                 4000+2000*(lambda==27), M, cv.plot=TRUE)
    title(title.string(M,lambda,lambda,seed))
    test.log.prob (tst.y, mlp.predict(tst.X,r))
  }
}

# Try some other possibilities.

seeds <- 1:2

lambda1 <- 9
lambda2 <- 9
M <- 8

for (seed in seeds)
{
  cat("\n")
  cat(title.string(M,lambda1,lambda2,seed),"\n")
  r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda1, lambda2, 1:nrow(val), 
               3000, M, cv.plot=TRUE)
  title(title.string(M,lambda1,lambda2,seed))
  test.log.prob (tst.y, mlp.predict(tst.X,r))
}

lambda1 <- 6
lambda2 <- 12
M <- 16

for (seed in seeds)
{
  cat("\n")
  cat(title.string(M,lambda1,lambda2,seed),"\n")
  r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda1, lambda2, 1:nrow(val), 
               4000, M, cv.plot=TRUE)
  title(title.string(M,lambda1,lambda2,seed))
  test.log.prob (tst.y, mlp.predict(tst.X,r))
}

lambda1 <- 12
lambda2 <- 6
M <- 16

for (seed in seeds)
{
  cat("\n")
  cat(title.string(M,lambda1,lambda2,seed),"\n")
  r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda1, lambda2, 1:nrow(val), 
               4000, M, cv.plot=TRUE)
  title(title.string(M,lambda1,lambda2,seed))
  test.log.prob (tst.y, mlp.predict(tst.X,r))
}

lambda1 <- 16
lambda2 <- 4
M <- 16

for (seed in seeds)
{
  cat("\n")
  cat(title.string(M,lambda1,lambda2,seed),"\n")
  r <- mlp.cv (trn.y, trn.X, 0.0001, 0.0002, lambda1, lambda2, 1:nrow(val), 
               4000, M, cv.plot=TRUE)
  title(title.string(M,lambda1,lambda2,seed))
  test.log.prob (tst.y, mlp.predict(tst.X,r))
}

# Close plot file.

dev.off()