Как построить границу решения классификатора k-ближайшего соседа из элементов статистического обучения?

31

Я хочу создать сюжет, описанный в книге ElemStatLearn «Элементы статистического обучения: сбор данных, вывод и прогноз. Второе издание» Тревора Хасти, Роберта Тибширани и Джерома Фридмана. Сюжет:

введите описание изображения здесь

Мне интересно, как я могу получить этот точный график R, особенно обратите внимание на графику сетки и расчеты, чтобы показать границу.

littleEinstein
источник
1
@StasK: да, это так. Как сформировать сюжет? Не могли бы вы помочь? Большое спасибо!
маленький Эйнштейн

Ответы:

35

Чтобы воспроизвести этот рисунок, вам необходимо установить пакет ElemStatLearn в вашей системе. Искусственный набор данных был сгенерирован mixture.example()как указано @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()

Все, кроме трех последних команд, взяты из интерактивной справки для mixture.example. Обратите внимание, что мы использовали тот факт, что expand.gridего выходные данные будут упорядочены путем изменения xвначале, что дополнительно позволяет индексировать (по столбцу) цвета в prob15матрице (размером 69x99), которая содержит долю голосов за класс победы для каждой координаты решетки ( px1, px2).

введите описание изображения здесь

хл
источник
+1. Благодарность! Мне также интересно, как генерировать данные, как описано в тексте «разоблачить оракула». Не могли бы вы также добавить это вместо использования данных с веб-сайта?
маленький Эйнштейн
@littleEinstein Вы имеете в виду то, что указано в онлайн-справке для mixture.example? Посмотрите на настройку симуляции ниже строки, начиная с # Reproducing figure 2.4, page 17 of the book:примера.
хл
пожалуйста, дайте мне знать ссылку? Я не могу это найти.
маленький Эйнштейн
Извините @littleEinstein, но есть кое-что, что я, вероятно, упускаю. Это всего лишь вопрос ввода help(mixture.example)или example(mixture.example)запроса R (после загрузки требуемого пакета с помощью library(ElemStatLearn)). Код для генерации искусственного набора данных (не для генерации Рис. 2.4) написан простым R в разделе Пример.
ЧЛ
1
Кстати, я только что наткнулся на блог @ Шейна, где он использовал ggplotдля подобных целей. Проверьте это: ESL 2.1: Линейная регрессия против KNN .
ЧЛ
7

Я самостоятельно изучаю ESL и пытаюсь проработать все примеры, представленные в книге. Я только что сделал это, и вы можете проверить код R ниже:

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 = rx1[1], to = rx1[2], by = 0.1 )
px2 <- seq(from = rx2[1], to = rx2[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()
Даоин Лин
источник
1
Чтобы ввести код здесь без этого, вы можете выделить текст, который является кодом, а затем нажать кнопку «Код» в верхней части страницы. Это в ряду иконок / кнопок. Код один выглядит как фигурные скобки.
Питер Флом - Восстановить Монику
Re: «Как вставить блок кода R». У вас есть доступ к небольшой строке меню при редактировании вашего сообщения.
ЧЛ
Кроме того, если вы не используете редактор, который может легко создавать отступы для блоков кода, я думаю, вы будете рады перейти на него. Например, в Rstudio выберите код и нажмите вкладку, чтобы сделать отступ, в VIM вы можете это сделать 5>>и т. Д.
Mark