KNN模型

根据KNN模型，使用R完成了KNN函数的实现代码，并与R中原先的KNN函数进行比较，最后把10折交叉检验的KNN模型与线性回归、二次回归、贝叶斯规则进行了比较

set.seed(012311)

Data Generation

Generate Centers

p = 2;      
csize = 10;     # number of centers
sigma = 1;      # sd for generating the centers 
m1 = matrix(rnorm(csize*p), csize, p)*sigma + 
  cbind( rep(1, csize), rep(0, csize))
m0 = matrix(rnorm(csize*p), csize, p)*sigma + 
  cbind( rep(0, csize), rep(1, csize))

Generate Data

sim_params = list(
 csize = 10,      # number of centers
 p = 2,           # dimension
 s = sqrt(1/5),   # standard deviation for generating data
 n = 100,         # training size per class
 N = 5000,        # test size per class
 m0 = m0,         # 10 centers for class 0
 m1 = m1         # 10 centers for class 1
)
generate_sim_data = function(sim_params){
  p = sim_params$p
  s = sim_params$s 
  n = sim_params$n 
  N = sim_params$N 
  m1 = sim_params$m1 
  m0 = sim_params$m0
  csize = sim_params$csize
  
  id1 = sample(1:csize, n, replace = TRUE);
  id0 = sample(1:csize, n, replace = TRUE);
  traindata = matrix(rnorm(2*n*p), 2*n, p)*s + rbind(m1[id1,], m0[id0,])
  Ytrain = factor(c(rep(1,n), rep(0,n)))
  shuffle_row_id = sample(1:n)
  id1 = sample(1:csize, N, replace=TRUE);
  id0 = sample(1:csize, N, replace=TRUE); 
  testdata = matrix(rnorm(2*N*p), 2*N, p)*s + rbind(m1[id1,], m0[id0,])
  Ytest = factor(c(rep(1,N), rep(0,N)))
  
  # Return the training/test data along with labels
  list(
  traindata = traindata,
  Ytrain = Ytrain,
  testdata = testdata,
  Ytest = Ytest
  )
}

Visualize Data

mydata = generate_sim_data(sim_params)
traindata = mydata$train
Ytrain = mydata$Ytrain
testdata = mydata$testdata
Ytest = mydata$Ytest
n = nrow(traindata)

mycol = rep("blue", n)
mycol[Ytrain==0] = "red"
plot(traindata[, 1], traindata[, 2], type = "n", xlab = "", ylab = "")
points(traindata[, 1], traindata[, 2], col = mycol);
points(m1[, 1], m1[, 2], pch = "+", cex = 2, col = "blue");    
points(m0[, 1], m0[, 2], pch = "+", cex = 2, col = "red");   
legend("bottomright", pch = c(1,1), col = c("blue", "red"), 
       legend = c("class 1", "class 0"))  

Part I

My own KNN coding

myknn <- function(train_set,train_Y,test_set,k) {
test_Y=vector()
for (i in c(1:dim(test_set)[1])) {
  dist_value=vector()
  for (j in c(1:dim(train_set)[1])){
    dist_value_matrix=dist(matrix(c(train_set[j,],test_set[i,]),nrow=2, ncol=dim(test_set)[2],byrow=TRUE),p=2)
    dist_value=c(dist_value,dist_value_matrix[1])
    arg_value=order(dist_value)[1:k]
    
  }
  Y=factor()
    for (i in c(1:k)){
      Y=c(Y,train_Y[arg_value[i]])
    }
    test_Y=c(test_Y,which.max(table(Y))-1)
  
  
} 
  return(test_Y)
}


bubbleSort<-function(vector) {
  n=length(vector)
  for (i in 1:(n-1)) {
    for (j in (i+1):n) {
      if(vector[i]>=vector[j]){
        temp = vector[i]
        vector[i] =vector[j]
        vector[j] = temp
        }
      }
    }
  return(vector)
}

Above is my own knn function.

About distance ties and voting ties

When it comes to voting ties and distance ties, I combine distance coefficient to the voting ties. I calculate the distance of all the data and sort it, then extract the first ‘k’ data which are closest to the predict point.

About the voting ties, the even k will not be allowed in the algorithm. But if k is even, in my code, it’s will be automatically labeled as 0.

I also have an idea about the voting ties, I can add a distance weight to the k data points, I can calculate the sum of all of the k data points with its weight. I think we can avoid the voting ties by this method.

Test

When k is 1.

Ytest_my=myknn(traindata,Ytrain,testdata,1)
table(Ytest,Ytest_my)

##      Ytest_my
## Ytest    0    1
##     0 4014  986
##     1  808 4192

library(class)
test.pred = knn(traindata, testdata, Ytrain, k = 1)
table(Ytest, test.pred)

##      test.pred
## Ytest    0    1
##     0 4014  986
##     1  808 4192

When k is 3.

Ytest_my=myknn(traindata,Ytrain,testdata,3)
table(Ytest,Ytest_my)

##      Ytest_my
## Ytest    0    1
##     0 4073  927
##     1  762 4238

library(class)
test.pred = knn(traindata, testdata, Ytrain, k = 3)
table(Ytest, test.pred)

##      test.pred
## Ytest    0    1
##     0 4073  927
##     1  762 4238

When k is 5.

Ytest_my=myknn(traindata,Ytrain,testdata,5)
table(Ytest,Ytest_my)

##      Ytest_my
## Ytest    0    1
##     0 4148  852
##     1  739 4261

library(class)
test.pred = knn(traindata, testdata, Ytrain, k = 5)
table(Ytest, test.pred)

##      test.pred
## Ytest    0    1
##     0 4148  852
##     1  739 4261

The above shows that my function ‘myknn’ matches the function ‘knn’.

Part II

Linear Regression

fit_lm_model = function(sim_data, verbose=FALSE) {
  
  # change Y from factor to numeric
  sim_data$Ytrain = as.numeric(sim_data$Ytrain) - 1
  sim_data$Ytest = as.numeric(sim_data$Ytest) - 1
  
  # fit a quadratic regression model
  model = lm(
    sim_data$Ytrain ~ 
      V1 + V2,
    as.data.frame(sim_data$traindata)
  )
  if (verbose) {
    print(summary(model))
  }
  
  decision_thresh = 0.5
  train_pred = as.numeric(model$fitted.values > decision_thresh)
 
  test_yhat = predict(
    model,
    newdata=as.data.frame(sim_data$testdata)
  )
  test_pred = as.numeric(test_yhat > decision_thresh)
  
  # return the mean classification errors on training/test sets
  list(
    train_error = sum(sim_data$Ytrain  != train_pred) / length(sim_data$Ytrain),
    test_error = sum(sim_data$Ytest  != test_pred) / 
      length(sim_data$Ytest)
  )
}

library(ggplot2)
lm_train=vector()
lm_test=vector()
for (i in c(1:50)){
  mydata = generate_sim_data(sim_params)
  lm_output = fit_lm_model(mydata)
  lm_train=c(lm_train,lm_output$train_error)
  lm_test=c(lm_test,lm_output$test_error)
}

Quadratic Regression

fit_qr_model_matrix = function(sim_data) {
  
  # change Y from factor to numeric
  sim_data$Ytrain = as.numeric(sim_data$Ytrain) - 1
  sim_data$Ytest = as.numeric(sim_data$Ytest) - 1
  
  train_matrix = cbind(sim_data$traindata, sim_data$traindata^2, sim_data$traindata[,1] * sim_data$traindata[,2])
  test_matrix = cbind(sim_data$testdata, sim_data$testdata^2, sim_data$testdata[,1] * sim_data$testdata[,2])
  
  # obtain quadratic regression coefs
  coefs = lm(sim_data$Ytrain ~ train_matrix)$coef
  train_yhat = coefs[1] + train_matrix %*% coefs[-1]
  test_yhat = coefs[1] + test_matrix %*% coefs[-1]
  
  decision_thresh = 0.5
  
  train_pred = as.numeric(train_yhat > decision_thresh)
  test_pred = as.numeric(test_yhat > decision_thresh)
  
  # return the mean classification errors on training/test sets
  list(
    train_error = sum(sim_data$Ytrain  != train_pred) / length(sim_data$Ytrain),
    test_error = sum(sim_data$Ytest  != test_pred) / 
      length(sim_data$Ytest)
  )
}

qr_train=vector()
qr_test=vector()
for (i in c(1:50)){
  mydata = generate_sim_data(sim_params)
  qr_output = fit_qr_model_matrix(mydata)
  qr_train=c(qr_train,qr_output$train_error)
  qr_test=c(qr_test,qr_output$test_error)
  
}

KNN classification with K chosen by 10-fold cross-validation

cvKNNAveErrorRate=function(K,traindata,Ytrain){
  
foldNum = 10
n = nrow(traindata)
foldSize = floor(n/foldNum)  
error = 0
myIndex = sample(1 : n)
for(runId in 1:foldNum){
  testSetIndex = ((runId-1)*foldSize + 1):(ifelse(runId == foldNum, n, runId*foldSize))
  testSetIndex = myIndex[testSetIndex]
  trainX = traindata[-testSetIndex, ]
  trainY = Ytrain[-testSetIndex]
  testX = traindata[testSetIndex, ]
  testY = Ytrain[testSetIndex]
  predictY = knn(trainX, testX, trainY, K)
  error = error + sum(predictY != testY) 
}
error = error / n

}

 cvKNN = function(traindata, Ytrain, foldNum) {
  n = nrow(traindata)
  foldSize = floor(n/foldNum)  
  KVector = seq(1, (nrow(traindata) - foldSize), 1)
  cvErrorRates = sapply(KVector, cvKNNAveErrorRate, traindata, Ytrain)
  result = list()
  result$bestK = max(KVector[cvErrorRates == min(cvErrorRates)])
  result$cvError = cvErrorRates[KVector == result$bestK]
  
  return(result)

}

cv_train=vector()
cv_test=vector()
for (i in c(1:50)){
  mydata = generate_sim_data(sim_params)
  cv_output = cvKNN(mydata$traindata,mydata$Ytrain,10)
  cv_train=c(cv_train,cv_output$cvError)
  k=cv_output$bestK
  test_pred=knn(mydata$traindata,mydata$testdata,mydata$Ytrain,k)
  test_error = sum(mydata$Ytest!= test_pred) / length(mydata$Ytest)
  cv_test=c(cv_test,test_error)
  
}

Bayes Rule

mixnorm = function(x, centers0, centers1, s){
  ## return the density ratio for a point x, where each 
  ## density is a mixture of normal with multiple components
  d1 = sum(exp(-apply((t(centers1) - x)^2, 2, sum) / (2 * s^2)))
  d0 = sum(exp(-apply((t(centers0) - x)^2, 2, sum) / (2 * s^2)))

  return (d1 / d0)
}

bayes_train_error=vector()
bayes_test_error=vector()
for (i in c(1:50)){
  mydata = generate_sim_data(sim_params)
  error1=0
  bayes_train_predict=vector()
  bayes_test_predict=vector()
  for (j in 1:200){
  a=mydata$traindata
  bayes_train=mixnorm(a[j,],m0,m1,sim_params$s)
  bayes_train_predict1=ifelse(bayes_train>=1,1,0)
  bayes_train_predict=c(bayes_train_predict,bayes_train_predict1)
  }
  error1 = sum(bayes_train_predict != mydata$Ytrain) 
  error1=error1/200
  bayes_train_error=c(bayes_train_error,error1)
  error2=0
  for (k in 1:10000){
  b=mydata$testdata
  bayes_test=mixnorm(b[k,],m0,m1,sim_params$s)
  bayes_test_predict1=ifelse(bayes_test>=1,1,0)
  bayes_test_predict=c(bayes_test_predict,bayes_test_predict1)
  }
  
  error2 = sum(bayes_test_predict != mydata$Ytest) 
  error2=error2/10000
  bayes_test_error=c(bayes_test_error,error2)
}

Plot

library(ggplot2)
par(pin = c(5,4))
boxplot(lm_train,lm_test,qr_train,qr_test,cv_train,cv_test,bayes_train_error,bayes_test_error,main="Errors for Four methods",col=terrain.colors(8))
#axis(2,seq(0,1,0.1),seq(0,1,0.1))
legend("topright", inset=0,   c("lm_train","lm_test","qr_train","qr_test","cv_train","cv_test","bayes_train","bayes_test"), fill=terrain.colors(8),bty='n')