source("ass2-funcs.r") pdf("ass2-plots.pdf",pointsize=6) # Read data and create vector of square-root responses and matrix with # original and computed covariates. d <- read.table("ass2-data.txt",head=T) n <- nrow(d) y <- sqrt(d$deaths) x <- cbind(d$temp, (1:n)/365.25, sin(2*pi*(1:n)/365.25), cos(2*pi*(1:n)/365.25)) # Define the two covariance functions that will be tried (differing in how # temperature is treated). covf1 <- function (x1, X2, h) { 10^2 + h[2]^2 * exp ( - ((x1[3]-X2[,3])^2 + (x1[4]-X2[,4])^2) / 2^2) + h[3]^2 * exp ( - abs(x1[2]-X2[,2]) / 0.1 ) + h[4]^2 * exp ( - abs(x1[2]-X2[,2]) / 0.02 ) + h[5]^2 * exp ( - (x1[1]-X2[,1])^2 / 20^2 ) } covf2 <- function (x1, X2, h) { 10^2 + h[2]^2 * exp ( - ((x1[3]-X2[,3])^2 + (x1[4]-X2[,4])^2) / 2^2) + h[3]^2 * exp ( - abs(x1[2]-X2[,2]) / 0.1 ) + h[4]^2 * exp ( - abs(x1[2]-X2[,2]) / 0.02 ) + h[5]^2 * exp ( - (x1[1]^3-X2[,1]^3)^2 / 20000^2 ) } # Use the first three years as training data. trn <- 1:(366+365+365) x_train <- x[trn,] y_train <- y[trn] # Estimate hyperparameters and make predictions with the two forms of the # covariance function. cat("With covf1, ") h1 <- gp_find_hypers (x_train, y_train, covf1, c (0.5, 0.1, 0.1, 0.1, 0.1)) cat("Hyperparameters with covf1:",round(h1,4),"\n\n") p1 <- gp_predict (x_train, y_train, x, covf1, h1) cat("With covf2, ") h2 <- gp_find_hypers (x_train, y_train, covf2, c (0.5, 0.1, 0.1, 0.1, 0.1)) cat("Hyperparameters with covf2:",round(h2,4),"\n\n") p2 <- gp_predict (x_train, y_train, x, covf2, h2) # Plot data and predictions with the two forms of the covariance function. par(mfrow=c(2,1)) plot(y,pch=20,col="red") points(p1,pch=20) abline(v=length(trn)+0.5) plot(y,pch=20,col="red") points(p2,pch=20) abline(v=length(trn)+0.5) # Compute average squared errors in predictions for last year, and for days # in the last year with temperature above 20 and below -5. cat("Average test error with covf1:",round(mean((y-p1)[-trn]^2),4),"\n") cat("Average test error with covf2:",round(mean((y-p2)[-trn]^2),4),"\n") cat("Average test error at temperature above 20 with covf1:", round(mean((y-p1)[-trn][x[-trn,1]>20]^2),4),"\n") cat("Average test error at temperature above 20 with covf2:", round(mean((y-p2)[-trn][x[-trn,1]>20]^2),4),"\n") cat("Average test error at temperature below -5 with covf1:", round(mean((y-p1)[-trn][x[-trn,1]< -5]^2),4),"\n") cat("Average test error at temperature below -5 with covf2:", round(mean((y-p2)[-trn][x[-trn,1]< -5]^2),4),"\n") # Make predictions with temperatures set one degree higher xx <- x xx[,1] <- x[,1] + 1 pp1 <- gp_predict (x_train, y_train, xx, covf1, h1) pp2 <- gp_predict (x_train, y_train, xx, covf2, h2) # Plot how change in prediction varies with temperature, for the # two forms of the covariance function, for training and test points. par(mfrow=c(2,2)) plot (x[trn,1], (pp1-p1)[trn], pch=20); abline(h=0) plot (x[trn,1], (pp2-p2)[trn], pch=20); abline(h=0) plot (x[-trn,1], (pp1-p1)[-trn], pch=20); abline(h=0) plot (x[-trn,1], (pp2-p2)[-trn], pch=20); abline(h=0) dev.off()