# (C) Petteri Pajunen 2006 # # # overfitting.R # # illustrate overfitting # iter <- 1000 X11(width=5,height=5) cat("Slope mu?") mu <- scan(n=1) beta <- 0 # intercept can be anything y0 <- rnorm(1,mean=beta,sd=1) y1 <- rnorm(1,mean=beta+mu,sd=1) plot(0,y0,lwd=3,type="p",col="green",xlim=c(-2,6),ylim=c(-5,10),xlab="x",ylab="y",main="GREEN: DATA, RED: LINEAR, BLUE: CONSTANT") points(0,beta,lwd=3) points(1,y1,lwd=3,col="green") points(1,beta+mu,lwd=3) xline <- -1:5 estmu <- y1-y0 estbeta <- y0 yline <- estmu*xline+estbeta lines(xline,matrix(1,1,7)*(0.5*(y0+y1)),col="blue") lines(xline,yline,col="red") scan(quiet=TRUE) linemse <- matrix(0,1,iter) linepred <- matrix(0,1,iter) constmse <- matrix(0,1,iter) constpred <- matrix(0,1,iter) for (i in 1:iter) { y0 <- rnorm(1,mean=beta,sd=1) y1 <- rnorm(1,mean=beta+mu,sd=1) linepred[i] <- 2*y1-y0 # predict at x=2 linemse[i] <- (2*y1-y0-beta-2*mu)^2 # square error constpred[i] <- 0.5*(y0+y1) constmse[i] <- (0.5*(y0+y1)-beta-2*mu)^2 ### square error } xp <- seq(min(min(linepred),min(constpred))-0.1,max(max(linepred),max(constpred))+0.5,by=0.5) linehist <- hist(linepred,breaks=xp,plot=FALSE) consthist <- hist(constpred,breaks=xp,plot=FALSE) plot(consthist$mids,consthist$intensities,type="l",col="blue",main="BLUE: CONSTANT, RED:LINEAR",xlab="PREDICTED VALUE AT x=2",ylab="DISTRIBUTION OF PREDICTIONS") lines(linehist$mids,linehist$intensities,col="red") points(beta+2*mu,0,lwd=5,col="black") text(beta+2*mu,max(linehist$intensities),labels=paste("LINE MSE: ",format(mean(linemse),digits=2)),col="red") text(beta+2*mu,max(consthist$intensities),labels=paste("CONST MSE: ",format(mean(constmse),digits=2)),col="blue") scan(quiet=TRUE) muvec <- seq(0.1,2,by=0.1) linemumse <- matrix(0,1,length(muvec)) constmumse <- matrix(0,1,length(muvec)) for (muiter in 1:length(muvec)) { mu <- muvec[muiter] cat("mu=",mu) cat("\n") for (i in 1:iter) { y0 <- rnorm(1,mean=beta,sd=1) y1 <- rnorm(1,mean=beta+mu,sd=1) linepred[i] <- 2*y1-y0 # predict at x=2 linemse[i] <- (2*y1-y0-beta-2*mu)^2 # square error constpred[i] <- 0.5*(y0+y1) constmse[i] <- (0.5*(y0+y1)-beta-2*mu)^2 ### square error } linemumse[muiter] <- mean(linemse) constmumse[muiter] <- mean(constmse) } plot(muvec,constmumse,lwd=2,col="blue",type="l",main="BLUE: CONSTANT, RED: LINEAR",xlab="TRUE VALUE OF mu",ylab="MSE PREDICTION ERROR") lines(muvec,linemumse,lwd=2,col="red") scan(quiet=TRUE) dev.off()