rm(list=ls())

library(class)
library(randomForest)

set.seed(456)

D = read.csv("dry-beans.csv",head=TRUE)     
C = D$Class
X = D[,1:16]     
n = ncol(X)
Y = rep(0,n)
Y[C == "DERMASON"] = 1
Y[C != "DERMASON"] = 0
n = length(Y)

I = sample(1:n,size=10000,replace=FALSE)
J = setdiff(1:n,I)
X = scale(X)
Xtrain = X[I,]
Xtest = X[J,]
Ytrain = Y[I]
Ytest = Y[J]

m = length(Ytest)
left = rep(0,20)
right = rep(0,20)
k = seq(1:20)
R = rep(0,20)
for(i in 1:20){
     out = knn(Xtrain,Xtest,Ytrain,k[i])
     z = out != Ytest
     R[i] = mean(z)
     se = sqrt(R[i]*(1-R[i])/m)
     left[i] = R[i] - 2*se
     right[i] = R[i] + 2*se
     }

LL=c(-3:-1,1:3)


pdf("knn_classification.pdf")
plot(1:20,R,pch=21,bg="red",ylim=c(0.03,.07),xlab="",ylab="",
             xaxt="n",cex.axis=1.3)
for(i in 1:20){
     segments(i,left[i],i,right[i],col="purple",lwd=4)
     }
points(1:20,R,pch=21,cex=1.2,bg="red")
LL = c(1,5,10,20)
values = c(1,5,10,20)
axis(1,LL,values, cex.axis=1.4,col.axis="black")
axis(1,at=14,"14",cex.axis=1.4,col.axis="red",col.lab="red")
mtext("Classification Error",side=2,line=2.3,cex=1.6)
mtext("Values of k",side=1,line=2.5,cex=1.6)
dev.off()


k = which.min(R)
Yhat = knn(Xtrain,Xtest,Ytrain,k)
print(table(Ytest,Yhat))

Ytrain = factor(Ytrain)
Ytest = factor(Ytest)
out = randomForest(Xtrain,Ytrain,xtest=Xtest,ytest=Ytest)
f = out$test$predicted
table(f,Ytest)

z= (Yhat != Ytest) - (f != Ytest)
Delta = mean(z)
s = sqrt(var(z))
left = Delta - 2*s/sqrt(m)
right = Delta + 2*s/sqrt(m)
cat("Delta = ",Delta,"\n")
print(c(left,right))