Como traçar o limite de decisão de um classificador de vizinhos k-mais próximos a partir de Elements of Statistical Learning?

31

Quero gerar o gráfico descrito no livro ElemStatLearn "Os elementos do aprendizado estatístico: mineração de dados, inferência e previsão. Segunda edição" de Trevor Hastie e Robert Tibshirani e Jerome Friedman. O enredo é:

insira a descrição da imagem aqui

Gostaria de saber como posso produzir esse gráfico exato R, principalmente os gráficos e cálculos da grade para mostrar o limite.

littleEinstein
fonte
1
@ Task: sim, é. Como gerar o enredo? Você poderia por favor ajudar? Muito Obrigado!
littleEinstein

Respostas:

35

Para reproduzir esta figura, você precisa ter o pacote ElemStatLearn instalado no seu sistema. O conjunto de dados artificial foi gerado com mixture.example()o apontado pelo @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()

Todos os comandos, com exceção dos três últimos, vêm da ajuda on-line mixture.example. Observe que usamos o fato de expand.gridorganizar sua saída variando xprimeiro, o que permite indexar (por coluna) as cores na prob15matriz (da dimensão 69x99), que mantém a proporção dos votos da classe vencedora para cada coordenada de treliça ( px1, px2).

insira a descrição da imagem aqui

chl
fonte
+1. obrigado! Também estou querendo saber como gerar os dados, conforme descrito no texto "expor o oráculo". Você também pode adicionar isso, em vez de usar os dados do site?
littleEinstein
@littleEinstein Você quer dizer o que é fornecido na ajuda on-line mixture.example? Veja a configuração da simulação abaixo da linha, começando # Reproducing figure 2.4, page 17 of the book:na seção de exemplo.
chl 24/01
você pode me informar o link? Eu não encontro.
littleEinstein
Desculpe @littleEinstein, mas provavelmente algo está faltando. É apenas uma questão de digitar help(mixture.example)ou example(mixture.example)no prompt do R (depois de carregar o pacote necessário library(ElemStatLearn)). O código para gerar o conjunto de dados artificial (não para gerar a Fig. 2.4) está escrito em R simples na seção Exemplo.
chl
1
Aliás, acabei de encontrar o blog do @ Shane, onde ele usava ggplotpara fins semelhantes. Confira: ESL 2.1: Regressão Linear vs. KNN .
chl
7

Sou autodidata em ESL e estou tentando trabalhar com todos os exemplos fornecidos no livro. Acabei de fazer isso e você pode conferir o código R abaixo:

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()
Daoying Lin
fonte
1
Para inserir o código aqui sem fazer isso, você pode realçar o texto que é o código e clicar no botão "código" na parte superior da página. Está em uma fileira de ícones / botões. O código parece um aparelho.
Peter Flom - Restabelece Monica
Re: "como colar um bloco de código R". Você tem acesso a uma barra de menu pequeno ao editar o seu post.
chl
Além disso, se você não estiver usando um editor que possa recuar blocos de código facilmente, acho que ficará feliz em mudar para um. Por exemplo, em rstudio selecionando o código e pressionar travessões guia-lo, no vim puder 5>>, etc.
Mark