# $Id: autompg.R,v 1.1 2007/10/02 10:40:14 kaip Exp $

#
# Copyright (c) 2007 Kai Puolamaki <Kai.Puolamaki@iki.fi>
#
# Permission to use, copy, modify, and distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#

AUTOS <- read.table("http://mlearn.ics.uci.edu/databases/auto-mpg/auto-mpg.data",
                 header=FALSE,na.strings="?")
colnames(AUTOS) <- c("mpg","cylinders","displacement","horsepower",
                     "weight","acceleration","model year","origin",
                     "car name")
AUTOS[,"origin"] <- factor(c("us","eu","jp")[AUTOS[,"origin"]])
AUTOS[,"mpg"] <- 235.4/AUTOS[,"mpg"]
colnames(AUTOS)[1] <- "lp100km"


AUTOS.lm <- lm(lp100km ~ horsepower,AUTOS)
AUTOS.lm10 <- lm(lp100km ~ horsepower+I(horsepower^2)+I(horsepower^3)+I(horsepower^4)+I(horsepower^5)+I(horsepower^6)+I(horsepower^7)+I(horsepower^8)+I(horsepower^9)+I(horsepower^10),AUTOS)

pdf("horsepower.pdf")
plot(AUTOS[,"horsepower"],AUTOS[,"lp100km"],
     pch=20,
     xlab="horsepower",ylab="fuel consumption (l/100km)",
     main="Linear regression")
a <- AUTOS.lm$coefficients["(Intercept)"]
b <- AUTOS.lm$coefficients["horsepower"]
abline(a=a,b=b,col="red")
text(50,25,sprintf("G(X)=%.2f+%.3f X",a,b),col="red",pos=4)
dev.off()

pdf("horsepower10.pdf")
plot(AUTOS[,"horsepower"],AUTOS[,"lp100km"],
     pch=20,
     xlab="horsepower",ylab="fuel consumption (l/100km)",
     main="Regression using degree 10 polynomial")
a <- AUTOS.lm10$coefficients["(Intercept)"]
b <- AUTOS.lm10$coefficients["horsepower"]
mm <- range(AUTOS[,"horsepower"],na.rm=TRUE)
xv <- mm[1]+(mm[2]-mm[1])*0:999/999
yv <- predict(AUTOS.lm10,data.frame(horsepower=xv))
lines(xv,yv,col="red")
text(50,25,sprintf("G(X)=%.1f%+.1f X+...",a,b),col="red",pos=4)
dev.off()