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

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.

-
Is it this one: www-stat.stanford.edu/~tibs/ElemStatLearn/datasets/… ? –  StasK Jan 24 '12 at 1:42

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). - +1. thanks! I am also wondering how to generate the data as described in the text "expose the oracle". Could you also please add that, instead of using the data from the website? – littleEinstein Jan 24 '12 at 19:32 @littleEinstein Do you mean what is given in the on-line help for mixture.example? Look at the simulation setup below line starting with # Reproducing figure 2.4, page 17 of the book: in the example section. – chl Jan 24 '12 at 20:12 can you please let me know the link? I cannot find it. – littleEinstein Jan 24 '12 at 20:51 Sorry @littleEinstein, but there's something I'm probably missing. It is just a matter of typing help(mixture.example) or example(mixture.example) at R prompt (after you load the required package with library(ElemStatLearn)). The code to generate the artificial dataset (not to generate Fig. 2.4) is written in plain R in the Example section. – chl Jan 24 '12 at 21:03 If I am not mistaken, the code there in the Example section actually references the data set, rather than generating it from scratch. I am talking about how to generate the dataset mixture.example itself. Really thank you for your help! – littleEinstein Jan 24 '12 at 23:31 show 2 more comments 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!!!

-
To enter code here without doing that, you can highlight the text that is code and then click on the "code" button near the top of the page. It's in a row of icons/buttons. The code one looks like braces. –  Peter Flom Oct 3 '12 at 10:17
Re: "how to paste a block of R code". You have access to a small menu bar when editing your post. –  chl Oct 3 '12 at 10:19

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.

-
The above code doesn't display any grid or boundary as requested, or did I miss something (I've tested the code, btw)? Anyway, shouldn't this be put as a new question? Cf. your last sentence. –  chl Sep 5 '12 at 21:32