#### Figure for Section 2, "Errors, In and Out of Sample" # Example: linear regression through the origin, with Gaussian noise # fix n, true slope, generate data, with X uniform on [0,1] n<-20; theta<-5 x<-runif(n); y<-x*theta+rnorm(n) # empirial risk = in-sample mean-squred error empirical.risk <- function(b) { mean((y-b*x)^2) } # generalization error of a line through the origin with slope b # EXERCISE: derive this formula true.risk <- function(b) { 1 + (theta-b)^2*(0.5^2+1/12) } # Plot the in-sample risk curve(Vectorize(empirical.risk)(x),from=0,to=2*theta, xlab="regression slope",ylab="MSE risk") # R trickery: the empirical.risk() function, as written, does not behalf # well when given a vector of slopes, and curve() wants its first argument to # be a function which can take a vector. Vectorize() turns its argument into # a function which can take a vector; writing the expression # Vectorize(empirical.risk)(x) # rather than just # Vectorize(empirical.risk) # helps curve() figure out where to pass its vector of points. curve(true.risk,add=TRUE,col="grey") # by contrast the true.risk() function works nicely with vectors. #### Figures for Section 3, "Over-Fitting and Model Selection" # Examples of training data # 20 standard-Gaussian X's x = rnorm(20) # Quadratic Y's y = 7*x^2 - 0.5*x + rnorm(20) # Initial plot of training data plus true regression curve plot(x,y) curve(7*x^2-0.5*x,col="grey",add=TRUE) # Fit polynomials and add them to the plot # Fit a constant (0th order polynomial), plot it y.0 = lm(y ~ 1) # We'd get the same constant by just doing mean(y), but fitting it as a # regression model means functions like residuals() and predict() are # available for use later, the same as our other models abline(h=y.0\$coefficients[1]) # Get evenly spaced points for pretty plotting of other models d = seq(min(x),max(x),length.out=200) # Fit polynomials of order 1 to 9 # It would be nicer if we let this run from 0 to 9, but R doesn't allow us # to do a polynomial of degree 0 for (degree in 1:9) { fm = lm(y ~ poly(x,degree)) # Store the results in models called y.1, y.2, through y.9 # The assign/paste trick here is often useful assign(paste("y",degree,sep="."), fm) # Plot them, with different line types lines(d, predict(fm,data.frame(x=d)),lty=(degree+1)) } # Calculate and plot in-sample errors mse.q = vector(length=10) for (degree in 0:9) { # The get() function is the inverse to assign() fm = get(paste("y",degree,sep=".")) mse.q[degree+1] = mean(residuals(fm)^2) } plot(0:9,mse.q,type="b",xlab="polynomial degree",ylab="mean squared error", log="y") # Plot the old curves with testing data x.new = rnorm(2e4) y.new = 7*x.new^2 - 0.5*x.new + rnorm(2e4) plot(x.new,y.new,xlab="x",ylab="y",pch=24,cex=0.1,col="blue") curve(7*x^2-0.5*x,col="grey",add=TRUE) abline(h=y.0\$coefficients[1]) d = seq(from=min(x.new),to=max(x.new),length.out=200) for (degree in 1:9) { fm = get(paste("y",degree,sep=".")) lines(d, predict(fm,data.frame(x=d)),lty=(degree+1)) } points(x,y) # Calculate and plot the out-of-sample errors gmse.q = vector(length=10) for (degree in 0:9) { # The get() function is the inverse to assign() fm = get(paste("y",degree,sep=".")) predictions = predict(fm,data.frame(x=x.new)) resids = y.new - predictions gmse.q[degree+1] = mean(resids^2) } plot(0:9,mse.q,type="b",xlab="polynomial degree", ylab="mean squared error",log="y",ylim=c(min(mse.q),max(gmse.q))) lines(0:9,gmse.q,lty=2,col="blue") points(0:9,gmse.q,pch=24,col="blue") ### Figures for section 4, "Cross-Validation" # Presumes the California housing data set from HW 1 has been loaded as a data # frame called calif # Display the first few lines of the first 3 columns head(calif[,1:3]) # Divide the whole data set in two at random which.half <- sample(rep(1:2,length.out=nrow(calif))) # Display the first few lines of the two halves head(calif[which.half==1,1:3]) head(calif[which.half==2,1:3]) # Fit the same simple regression to each half small.1 <- lm(MedianHouseValue~MedianIncome,data=calif[which.half==1,]) small.2 <- lm(MedianHouseValue~MedianIncome,data=calif[which.half==2,]) # Fit the same multiple regression to each half large.1 <- lm(MedianHouseValue~MedianIncome+MedianHouseAge, data=calif[which.half==1,]) large.2 <- lm(MedianHouseValue~MedianIncome+MedianHouseAge, data=calif[which.half==2,]) # Function to calculate in-sample mean-squared error in.sample.mse <- function(fit) { mean(fit\$residuals^2) } in.sample.mse(small.1) in.sample.mse(large.1) in.sample.mse(small.2) in.sample.mse(large.2) # Function to calculate the error of a given function on a given half of # the data # EXERCISE: does cross.sample.mse(small.1,1) == in.sample.mse(small.1)? # EXERCISE: Should they be equal? cross.sample.mse <- function(fit,half) { test <- calif[which.half != half,] return(mean((test\$MedianHouseValue - predict(fit,newdata=test))^2)) } # Find errors of each model on the other half of the data cross.sample.mse(small.1,2) cross.sample.mse(small.2,1) cross.sample.mse(large.1,2) cross.sample.mse(large.2,1) # Extracting the coefficients and averaging the cross-prediction scores # is left as an easy EXERCISE