# STA 414/2104, Spring 2006, Assignment #4.
#
# R. M. Neal
#
# Uses the Gaussian process regression functions from gp.r.

source("gp.r")

# The additive covariance function.

covf1 <- function (x1, x2, h)
{
  eta <- h[2]
  rho <- h[3:length(h)]

  10^2 + eta^2 * sum (exp(-((x1-x2)/rho)^2))
}

# The non-additive covariance function.

covf2 <- function (x1, x2, h)
{
  eta <- h[2]
  rho <- h[3:length(h)]

  10^2 + eta^2 * exp ( - sum( ((x1-x2)/rho)^2 ) )
}


# Read the data for assignment 4.

x.train4 <- as.matrix(read.table("a4-x-train",head=F))
y.train4 <- scan("a4-y-train")
x.test4 <- as.matrix(read.table("a4-x-test",head=F))
y.test4 <- scan("a4-y-test")
x.grid4 <- as.matrix(read.table("a4-x-grid",head=F))

# Make predictions with the additive covariance function.

h4.1 <- gp.find.hypers (x.train4, y.train4, covf1, c(1,1,1,1))
p4.1 <- gp.predict (x.train4, y.train4, x.test4, covf1, h4.1)
g4.1 <- gp.predict (x.train4, y.train4, x.grid4, covf1, h4.1)

# Make predictions with the non-additive covariance function.

h4.2 <- gp.find.hypers (x.train4, y.train4, covf2, c(1,1,1,1))
p4.2 <- gp.predict (x.train4, y.train4, x.test4, covf2, h4.2)
g4.2 <- gp.predict (x.train4, y.train4, x.grid4, covf2, h4.2)

# Report average squared errors.

cat ("A4 additive average squared error: ",mean((y.test4-p4.1)^2),"\n")
cat ("A4 non-additive average squared error: ",mean((y.test4-p4.2)^2),"\n")

# Report log probability of data.

cat ("A4 additive log probability: ",
       gp.log.prob (x.train4, y.train4, covf1, h4.1), "\n")
cat ("A4 non-additive log probability: ",
       gp.log.prob (x.train4, y.train4, covf2, h4.2), "\n")

# Read the data for assignment 3.

x.train3 <- as.matrix(read.table("a3-x-train",head=F))
y.train3 <- scan("a3-y-train")
x.test3 <- as.matrix(read.table("a3-x-test",head=F))
y.test3 <- scan("a3-y-test")
x.grid3 <- as.matrix(read.table("a3-x-grid",head=F))

# Make predictions with the additive covariance function.

h3.1 <- gp.find.hypers (x.train3, y.train3, covf1, c(1,1,1,1))
p3.1 <- gp.predict (x.train3, y.train3, x.test3, covf1, h3.1)
g3.1 <- gp.predict (x.train3, y.train3, x.grid3, covf1, h3.1)

# Make predictions with the non-additive covariance function.

h3.2 <- gp.find.hypers (x.train3, y.train3, covf2, c(1,1,1,1))
p3.2 <- gp.predict (x.train3, y.train3, x.test3, covf2, h3.2)
g3.2 <- gp.predict (x.train3, y.train3, x.grid3, covf2, h3.2)

# Report average squared errors.

cat ("A3 additive average squared error: ",mean((y.test3-p3.1)^2),"\n")
cat ("A3 non-additive average squared error: ",mean((y.test3-p3.2)^2),"\n")

# Report log probability of data.

cat ("A3 additive log probability: ",
       gp.log.prob (x.train3, y.train3, covf1, h3.1), "\n")
cat ("A3 non-additive log probability: ",
       gp.log.prob (x.train3, y.train3, covf2, h3.2), "\n")


# Produce contour plots.

postscript("ass4-plt.ps",pointsize=10,horiz=F,height=8,width=7)
par(mfrow=c(2,2))

contour (seq(0,1,0.02), seq(0,1,0.02), matrix(g4.1,51,51), 
         xlab="A4 data, additive covariance")

contour (seq(0,1,0.02), seq(0,1,0.02), matrix(g4.2,51,51), 
         xlab="A4 data, non-additive covariance")

contour (seq(0,1,0.02), seq(0,1,0.02), matrix(g3.1,51,51), 
         xlab="A3 data, additive covariance")

contour (seq(0,1,0.02), seq(0,1,0.02), matrix(g3.2,51,51), 
         xlab="A3 data, non-additive covariance")

dev.off()