How to plot decision boundary of a k-nearest neighbor classifier from Elements of Statistical Learning?

Question

I want to generate the plot described in the book ElemStatLearn "The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Second Edition" by Trevor Hastie & Robert Tibshirani& Jerome Friedman. The plot is:

I am wondering how I can produce this exact graph in R, particularly note the grid graphics and calculation to show the boundary.

Answer 1

To reproduce this figure, you need to have the ElemStatLearn package installed on you system. The artificial dataset was generated with mixture.example() as pointed out by @StasK.

library(ElemStatLearn)
require(class)
x <- mixture.example$x
g <- mixture.example$y
xnew <- mixture.example$xnew
mod15 <- knn(x, xnew, g, k=15, prob=TRUE)
prob <- attr(mod15, "prob")
prob <- ifelse(mod15=="1", prob, 1-prob)
px1 <- mixture.example$px1
px2 <- mixture.example$px2
prob15 <- matrix(prob, length(px1), length(px2))
par(mar=rep(2,4))
contour(px1, px2, prob15, levels=0.5, labels="", xlab="", ylab="", main=
        "15-nearest neighbour", axes=FALSE)
points(x, col=ifelse(g==1, "coral", "cornflowerblue"))
gd <- expand.grid(x=px1, y=px2)
points(gd, pch=".", cex=1.2, col=ifelse(prob15>0.5, "coral", "cornflowerblue"))
box()

All but the last three commands come from the on-line help for mixture.example. Note that we used the fact that expand.grid will arrange its output by varying x first, which further allows to index (by column) colors in the prob15 matrix (of dimension 69x99), which holds the proportion of the votes for the winning class for each lattice coordinates (px1,px2).

Answer 2

I'm self-learning ESL and trying to work through all examples provided in the book. I just did this and you can check the R code below:

library(MASS)
# set the seed to reproduce data generation in the future
seed <- 123456
set.seed(seed)

# generate two classes means
Sigma <- matrix(c(1,0,0,1),nrow = 2, ncol = 2)
means_1 <- mvrnorm(n = 10, mu = c(1,0), Sigma)
means_2 <- mvrnorm(n = 10, mu = c(0,1), Sigma)

# pick an m_k at random with probability 1/10
# function to generate observations
genObs <- function(classMean, classSigma, size, ...)
{
  # check input
  if(!is.matrix(classMean)) stop("classMean should be a matrix")
  nc <- ncol(classMean)
  nr <- nrow(classMean)
  if(nc != 2) stop("classMean should be a matrix with 2 columns")
  if(ncol(classSigma) != 2) stop("the dimension of classSigma is wrong")

  # mean for each obs
    # pick an m_k at random
  meanObs <- classMean[sample(1:nr, size = size, replace = TRUE),]
  obs <- t(apply(meanObs, 1, function(x) mvrnorm(n = 1, mu = x, Sigma = classSigma )) )
  colnames(obs) <- c('x1','x2')
  return(obs)
}


obs100_1 <- genObs(classMean = means_1, classSigma = Sigma/5, size = 100)
obs100_2 <- genObs(classMean = means_2, classSigma = Sigma/5, size = 100)

# generate label
y <- rep(c(0,1), each = 100)

# training data matrix
trainMat <- as.data.frame(cbind(y, rbind(obs100_1, obs100_2)))

# plot them
library(lattice)
with(trainMat, xyplot(x2 ~ x1,groups = y, col=c('blue', 'orange')))

# now fit two models

# model 1: linear regression
lmfits <- lm(y ~ x1 + x2 , data = trainMat)

# get the slope and intercept for the decision boundary
intercept <- -(lmfits$coef[1] - 0.5) / lmfits$coef[3]
slope <- - lmfits$coef[2] / lmfits$coef[3]

# Figure 2.1
xyplot(x2 ~ x1, groups = y, col = c('blue', 'orange'), data = trainMat,
       panel = function(...)
       {
        panel.xyplot(...)
        panel.abline(intercept, slope)
        },
       main = 'Linear Regression of 0/1 Response')    

# model2: k nearest-neighbor methods
library(class)
# get the range of x1 and x2
rx1 <- range(trainMat$x1)
    rx2 <- range(trainMat$x2)
# get lattice points in predictor space
px1 <- seq(from = -1.6, to = 4.0, by = 0.1 )
px2 <- seq(from = -2, to = 3.2, by = 0.1 )
xnew <- expand.grid(x1 = px1, x2 = px2)

# get the contour map
knn15 <- knn(train = trainMat[,2:3], test = xnew, cl = trainMat[,1], k = 15, prob = TRUE)
prob <- attr(knn15, "prob")
prob <- ifelse(knn15=="1", prob, 1-prob)
prob15 <- matrix(prob, nrow = length(px1), ncol = length(px2))

# Figure 2.2
par(mar = rep(2,4))
contour(px1, px2, prob15, levels=0.5, labels="", xlab="", ylab="", main=
    "15-nearest neighbour", axes=FALSE)
points(trainMat[,2:3], col=ifelse(trainMat[,1]==1, "coral", "cornflowerblue"))
points(xnew, pch=".", cex=1.2, col=ifelse(prob15>0.5, "coral", "cornflowerblue"))
box()

btw, could anyone tell me how to paste a block of R code to here without entering four white spaces for each line every time??? it's kinda pain to do this!!!

Answer 3

here is my code to generate the data and graphs. It's different from ElemStatLearn.

library(MASS)
n = 10
mu1 = c(1,0)
Sigma1 = matrix(c(1,0,0,1), c(2,2))
center1 = mvrnorm(n, mu = mu1, Sigma = Sigma1)
mu2 = c(0,1)
Sigma2 = Sigma1
center2 = mvrnorm(n, mu = mu2, Sigma = Sigma2)
mu1.series = ceiling(runif(n^2, max =10))
mu2.series = ceiling(runif(n^2, max =10))
Sigma.pt = Sigma1/5

blue.pt = t(apply(center1[mu1.series, ], 1, function(m){mvrnorm(1, m, Sigma.pt)}))
orange.pt = t(apply(center2[mu2.series, ], 1, function(m){mvrnorm(1, m, Sigma.pt)}))

x.min = min(c(blue.pt[, 1], orange.pt[, 1]))
x.max = max(c(blue.pt[, 1], orange.pt[, 1]))
y.min = min(c(blue.pt[, 2], orange.pt[, 2]))
y.max = max(c(blue.pt[, 1], orange.pt[, 2]))
x.range = seq(from = x.min, to = x.max, by = 0.1)
y.range = seq(from = y.min, to = y.max, by = 0.1)
x.length = length(x.range)
y.length = length(y.range)
pt.all = cbind(rep(x.range, y.length), rep(y.range, x.length))
thre = 0.2


plot(c(blue.pt[, 1], orange.pt[, 1]), c(blue.pt[, 2], orange.pt[, 2]), type = "n",
 main = "Bayes Optimal Classifer", xlab = "X-axis", ylab = "Y-axis")
points(blue.pt, col = 'blue')
points(orange.pt, col = 'orange')

May I ask do you know how to calculate the Bayes decision boundary in Figure 2.5? Thank you very much.