# code based on the files "Code regarding section 3.2 (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa)" # and "Code regarding section 3.3-3.5" (copyright 2003, 2004, 2012 A.Azzalini and B.Scarpa) # read data dataset <- read.table("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", header = TRUE) attach(dataset) n <- nrow(dataset) y.range <- range(c(y.yesterday, y.tomorrow)) # define f(x) as polynomials polinomio <- function(p) if(p == 0) lm(y.yesterday ~ 1) else lm(y.yesterday ~ poly(x, p)) polinomi <- apply(cbind(0:23), 1, polinomio) # grid to plot curves x0 <- data.frame(x = sort(unique(c(x, seq(min(x), max(x), length = 200))))) ### training set ### # figure 3.1 plot(x, y.yesterday, ylim = y.range, pch = 16, ylab = "y", cex = 2) points(x, y.tomorrow, ylim = y.range, pch = 16, ylab = "y", cex = 2) # figure 3.2 # a plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 3rd degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[4]], newdata = x0), lwd = 3) # b plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 6th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[7]], newdata = x0), lwd = 3) # c plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 12th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[13]], newdata = x0), lwd = 3) # d plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 18th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[19]], newdata = x0), lwd = 3) # e plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 24th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[24]], newdata = x0), lwd = 3) # f # plot(x, y.yesterday ,ylim = y.range, col = 2, pch = 16, # main = "Data and (n-1)th Degree Polynomial") # lines(x0$x, predict(polinomi[[n]], newdata = x0)) # x.est <- sapply(1:n, function(x, data) data^x, data = x) # it does not work, look at the picture in the book instead (taken with an older version of R) # all in one par(mfrow = c(2, 2)) plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 3rd degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[4]], newdata = x0), lwd = 3) plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 6th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[7]], newdata = x0), lwd = 3) plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 12th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[13]], newdata = x0), lwd = 3) plot(x, y.yesterday, ylim = y.range, col = 2, pch = 16, xlab = "", main = "data and 24th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[24]], newdata = x0), lwd = 3) par(mfrow = c(1, 1)) par(mfrow = c(1, 1)) # figure 3.3a dev <- as.numeric(lapply(polinomi, deviance)) plot(0:(23), dev, type = "b", xlab = "p", ylab = "Deviance", ylim = c(0,max(dev)), cex = 2, lwd = 2) # figure 3.3b R2 <- function(obj) summary(obj)$r.squared r2 <- as.numeric(lapply( polinomi, R2)) plot(0:(23), r2, type="b", xlab="p", ylab=expression(R^2), cex = 2, lwd = 2) ### test set ### plot(x, y.yesterday, ylim = y.range, pch = 16, ylab = "y", cex = 2, col = 2) points(x, y.tomorrow, ylim = y.range, pch = 17, ylab = "y", cex = 2, col = 4) # figure 3.4a plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 3rd degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[4]], newdata = x0), lwd = 3) # figure 3.4b plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 6th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[7]], newdata = x0), lwd = 3) # figure 3.4c plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 12th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[13]], newdata = x0), lwd = 3) # figure 3.4d plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 18th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[19]], newdata=x0), lwd = 3) # figure 3.4e plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 24th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[25]], newdata = x0), lwd = 3) # figure 3.4f # plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", # main="data and (n-1)th degree polynomial", font.main=1, cex = 2) # lines(x0$x, predict(polinomi[[n]], newdata=x0), lwd = 3) # it does not work, look at the picture in the book instead (taken with an older version of R) # all in one par(mfrow = c(2, 2)) plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 3rd degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[4]], newdata = x0), lwd = 3) plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 6th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[7]], newdata = x0), lwd = 3) plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 12th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[13]], newdata = x0), lwd = 3) plot(x, y.tomorrow, ylim = y.range, col = 4, pch = 17, xlab = "", main = "data and 24th degree polynomial", font.main = 1, cex = 2) lines(x0$x, predict(polinomi[[25]], newdata = x0), lwd = 3) par(mfrow = c(1, 1)) # figure 3.5a y.hat <- lapply( polinomi, predict) devianza <- function( y.hat) sum((y.tomorrow - y.hat)^2) dev2 <- as.numeric(lapply(y.hat, devianza)) plot(0:23, dev2, type = "b", xlab = "p", ylab = "deviance", ylim=c(0,max(dev)), cex = 2, lwd = 2) # figure 3.5b r2 <- 1 - dev2 / dev[1] plot(0:23, r2, type = "b", xlab = "p", ylab = expression(R^2), ylim = c(0, 1), cex = 2, lwd = 2) polynomial <- function(p) lm(y.yesterday ~ poly(x,p)) #it only works until polynomial of degree 23 (up to 29) polynomials <- apply(cbind(1:23), 1, polynomial) error.train <- sapply(polynomials, deviance) rss <- function(p) sum((y.tomorrow - predict(polynomials[[p]]))^2) error.test <- apply(cbind(1:23), 1, rss) out <- cbind(1:23, error.train, error.test) colnames(out) <- c('complexity', 'train.error', 'test.error') out # bias-variance trade-off rss60<- function(p) sum((y.tomorrow - predict(polynomials[[p]]))^2)+ sum((y.yesterday - predict(polynomials[[p]]))^2) Error <- apply(cbind(1:23), 1, rss) source("f_true.R") polo.f.true <- function(p) lm(f.true ~ poly(x,p)) poli.f.true <- apply(cbind(1:23), 1, polo.f.true) bias2 <- function(p) mean((fitted(poli.f.true[[p]]) - f.true)^2) bias2.oss <- apply(cbind(1:23), 1, bias2) var.oss <- sqm.true^2 * (1:23)/30 err <- bias2.oss + var.oss # figure 3.6 plot(1:23, err, type = "b", cex = 2, pch = 16, xlab = "Model complexity", ylab = "Error") # figure 3.7 plot(bias2.oss, type = "b", col = 4, ylab = "Error", xlab = "Model complexity", cex = 2, pch = 16) lines(var.oss, type = "b", col = 2, cex = 2, pch = 17) lines(err, type = "b", cex = 2, pch = 18) legend('topright', legend = c("Bias^2", "Variance", "Total"), col = c(4, 2, 1), pch = c(16, 17, 18)) # figure 3.8 plot(error.train, type = "b", col = 2, xlab = "Model complexity", ylab = "Deviance", cex = 2, pch = 16) lines(error.test, type = "b",col = 4, cex = 2, pch = 17) text(x = c(19, 19), y = c(0.0015, 0.0065), c("Yesterday data", "Tomorrow data")) ### cross-validation ### # it only works until polynomial of degree 23 (up to 29) # figure 3.9 y <- c(y.yesterday, y.tomorrow) xx <- c(x, x) RSS <- numeric(23) for (p in 0:22) { m <- lm(y ~ poly(xx, p+1), x = TRUE) X <- m$x P <- X %*% solve(t(X)%*% X) %*% t(X) # projection matrix p.ii <- diag(P) RSS.pp <- sum(residuals(m)^2 / (1 - p.ii)^2) RSS.p <- 0 RSS[p + 1] <- RSS.pp } plot(0:22, RSS, type = "b", xlab = "Model complexity", pch = 16, ylab = "Error", cex = 2, lwd = 3)