rm(list=ls())


library(randomForest)

D = read.table("Auto.data",head=TRUE)
D = D[,-c(8,9)]
D[,4] = as.numeric(D[,4])
I = complete.cases(D)
D = D[I,]

X = D[,2:7]
Y = D[,1]
X = as.matrix(X)

n = nrow(D)
m = round(4*n/5)
I = sample(1:n,size=m,replace=FALSE)

Xtrain = X[I,]
Ytrain = Y[I]
Xtest = X[-I,]
Ytest = Y[-I]

## linear model
train = data.frame(X=Xtrain,Y=Ytrain)
out = lm(Y ~ .,data=train)
test = data.frame(X=Xtest)
p = predict(out,newdata=test)
error = mean((Ytest-p)^2)
plot(Ytest,p,xlim=c(5,50),ylim=c(5,50))
print(error)
e1 = (Ytest-p)^2

## forest
out = randomForest(x=Xtrain,y=Ytrain,xtest=Xtest,ytest=Ytest,
                   importance=TRUE,keep.forest=TRUE)
p = predict(out,newdata=Xtest)
names(p) = NULL
points(Ytest,p,col="red")
abline(a=0,b=1)
error = mean((Ytest-p)^2)
print(error)
e2 = (Ytest-p)^2

e = e1-e2
d = mean(e)
s = sqrt(var(e))
m = length(e)
print(c(d - 2*s/sqrt(m),d + 2*s/sqrt(m)))