Commit 5f316d80 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files

Adiciona mais conteúdo para aula de SVM.

parent c2a6a37a
......@@ -31,6 +31,8 @@ navbar:
href: tutorials/03-regularization.html
- text: "Métodos de árvores"
href: tutorials/04-regression-trees.html
- text: "Máquina de vetores de suporte"
href: tutorials/05-support-vector-machine.html
- text: "Scripts"
icon: fa-file-text
href: scripts/
......
---
title: "Máquinas de vetores de suporte"
author: Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira
date: 2017-10-10
#bibliography: ../config/Refs.bib
#csl: ../config/ABNT-UFPR-2011-Mendeley.csl
---
```{r, include = FALSE}
source("../config/setup.R")
opts_chunk$set(
cache = FALSE,
message = FALSE,
warning = FALSE)
```
# Classificador de vetores de suporte
```{r, message = FALSE}
rm(list = objects())
library(lattice)
library(latticeExtra)
```
## Simulando um caso linearmente separável
```{r}
#-----------------------------------------------------------------------
# Gerando um conjunto de dados artifical.
s <- seq(0, 1, length.out = 7)
da <- expand.grid(x1 = s, x2 = s)
# Se acima da linha 1:1 + 0.5, então 1, senão -1.
da$y <- with(da, ifelse(0.1 + x1 - x2 < 0, "blue", "red"))
with(da, plot(x2 ~ x1, col = y, asp = 1))
abline(a = 0.1, b = 1, col = "magenta")
```
Uma solução construída usando o MATLAB:
<https://oceandatamining.sciencesconf.org/conference/oceandatamining/TP_Linear_SVM.pdf>.
## Transformando em não linearmente separável
```{r}
#-----------------------------------------------------------------------
# Versão com um jitter nos pontos para dar uma misturada.
set.seed(1234)
a <- 0.1
db <- transform(da,
x1 = x1 + runif(nrow(da), -a, a),
x2 = x2 + runif(nrow(da), -a, a))
with(db, plot(x2 ~ x1, col = y, asp = 1))
abline(a = 0.1, b = 1, col = "magenta")
```
## Otimização quadrática com restrições
Veja aqui
<http://members.cbio.mines-paristech.fr/~thocking/mines-course/2011-04-01-svm/svm-qp.pdf>
uma tradução do problema de classificador de vetores de suporte usando
otimização convexa com o pacote
[`quadprog`](https://cran.r-project.org/web/packages/quadprog).
A implementação de classificador de vetores de suporte com o `quadprog`
foi apresentada no stackoverflow:
<https://stackoverflow.com/questions/33088363/how-to-use-r-package-quadprog-to-solve-svm>.
```{r}
# Pacote para problemas de otimização convexa com restrições lineares.
library(quadprog)
# Acesso à documentação.
# help(package = "quadprog", h = "html")
# Aplicando para os dados de duas especies de iris com as variáveis de
# comprimento apenas.
names(iris)
da <- droplevels(subset(iris,
# Species != "versicolor",
# Species != "virginica",
Species != "setosa",
select = c(1, 3, 5)))
# Codifica categorias com -1 e 1.
y <- 2 * as.integer(da$Species) - 3
# Matriz de variáveis preditoras escalonadas.
X <- scale(as.matrix(da[, 1:2]))
# Cabeça e cauda das matrizes.
head(cbind(X, y))
tail(cbind(X, y))
# Diagrama de dispersão com indicação das classes.
plot(X, col = y + 3)
#--------------------------------------------
# Desenvolvendo a solução com o quadprog.
# Número de observações.
N <- nrow(X)
# Número de dimensões.
n_d <- ncol(X)
# Valor do parâmetro de regularização.
C <- 1
# Vetor que multiplica o vetor de parâmetros
# b = [\beta_0, \beta, \zeta].
d <- c(0,
rep(0, n_d),
rep(-C, N))
# Valor positivo pequeno mas não nulo para fazer a matriz D ser positiva
# definida (exigência da quadprog()).
eps <- 1e-10
D <- diag(c(eps,
rep(1, n_d),
rep(eps, N)))
# Matriz diagonal de dimensão N.
I_N <- diag(N)
# Concatena matrizes de zeros (dos betas) com a diagonal (dos epsilons).
# Representar a restrição zeta_i > 0, para todo i = 1, ..., N.
A_1 <- cbind(0,
matrix(0, ncol = n_d, nrow = N),
I_N)
dim(A_1)
# Representa a condição zeta_i + beta_0 * y_i + beta^T * x_i * y_i >= 1,
# para todo i = 1, ..., N.
A_2 <- as.matrix(cbind(as.matrix(y),
X * as.matrix(y)[, rep(1, n_d)],
I_N))
dim(A_2)
# Elimina nomes de linhas e colunas.
rownames(A_1) <- NULL
colnames(A_1) <- NULL
rownames(A_2) <- NULL
colnames(A_2) <- NULL
# Concatena as matrizes e vetores.
A <- t(rbind(A_1, A_2))
dim(A)
# image(A, asp = 1)
b_0 <- c(rep(0, N),
rep(1, N))
length(b_0)
# Passando os componentes para o otimizador.
results <- solve.QP(D, d, A, b_0)
# Valor da função objetivo na convergência.
results$value
# Extrai os parâmetros otimizados, i.e., [beta_0, beta, zeta].
b_optim <- results$solution
b_optim
# Hiperplano de separação.
beta_0 <- b_optim[1]
beta_0
beta <- b_optim[1 + (1:n_d)]
beta
# 1/sum(beta^2)
# Coeficientes.
zeta <- b_optim[(n_d + 1) + (1:N)]
zeta
# Vetores de suporte.
sp <- rbind(X[which(zeta > 0), ])
sp
nrow(sp)
# Predição a partir do modelo ajustado.
y_pred <- sign(apply(X,
MARGIN = 1,
FUN = function(x) {
beta_0 + sum(beta * as.vector(x))
}))
table(y, y_pred)
# Classificações incorretas.
i <- which(y * y_pred < 0)
length(i)
# Fazendo um grid para predição.
grid <- with(da,
expand.grid(
seq(min(X[, 1]), max(X[, 1]), length.out = 41),
seq(min(X[, 2]), max(X[, 2]), length.out = 41),
KEEP.OUT.ATTRS = FALSE))
names(grid) <- names(da)[1:2]
grid$y_pred <- sign(apply(as.matrix(grid),
MARGIN = 1,
FUN = function(x) {
beta_0 + sum(beta * as.vector(x))
}))
# Gráfico com pontos, classficações, fronteira e vetores de suporte.
plot(X,
col = y + 3,
asp = 1,
cex = 1.2,
pch = 19)
points(grid[, 1:2],
col = grid$y_pred + 3,
pch = 3)
abline(a = -beta_0/beta[2],
b = -beta[-2]/beta[2],
lty = 2,
col = "orange")
arrows(x0 = 0,
y0 = 0,
x1 = sp[, 1],
y1 = sp[, 2],
pch = 4,
cex = 2,
col = "magenta",
length = 0.05)
if (length(i)) {
points(rbind(X[i, ]),
pch = 19,
cex = 0.8,
col = "yellow")
}
```
# Uso de funções kernel
A ser feito.
# Implementações de máquinas de vetores de suporte
Uma revisão das implementações de máquina de vetores de suporte para o
ambiente R está disponível em
<https://www.jstatsoft.org/article/view/v015i09/v15i09.pdf>. Essa
revisão é de 2006, por isso, certamente foi feito
aperfeiçoamento/extensão dos recursos disponíveis além do surgimento de
novas implementações.
## Pacote `e1071`
Para fins didáticos, o mesmo problema resolvido com a implementação
baseada no pacote `quadprog` é revisitado usando implementações de SVM.
O pacote [`e1071`](https://cran.r-project.org/web/packages/e1071/)
possui a função `svm()` que dá acesso a implementação disponível na
[`libsvm`](https://www.csie.ntu.edu.tw/~cjlin/libsvm/) desenvolvida em
C++.
```{r}
library(e1071)
# Acessa a documentação da função.
# help(svm, help_type = "html")
m0 <- svm(x = X,
y = y,
scale = FALSE,
kernel = "linear",
type = "C-classification",
cost = 1)
m0
# Classe e métodos.
class(m0)
methods(class = class(m0))
# Tabela de confusão.
table(y, predict(m0))
# Classificações incorretas.
i <- which(y * as.integer(as.character(predict(m0))) < 0)
length(i)
# Vetores de suporte.
sp <- m0$SV
sp
# Fazendo a predição.
grid$y_svm <- predict(m0, newdata = grid[, 1:2])
# Gráfico com pontos, classficações, fronteira e vetores de suporte.
plot(X,
col = y + 3,
asp = 1,
cex = 1.2,
pch = 19)
points(grid[, 1:2],
col = 2 * as.integer(grid$y_svm),
pch = 3)
abline(a = -beta_0/beta[2],
b = -beta[-2]/beta[2],
lty = 2,
col = "orange")
arrows(x0 = 0,
y0 = 0,
x1 = sp[, 1],
y1 = sp[, 2],
pch = 4,
cex = 2,
col = "magenta",
length = 0.05)
if (length(i)) {
points(rbind(X[i, ]),
pch = 19,
cex = 0.8,
col = "yellow")
}
```
## Pacote `kernlab`
```{r}
library(kernlab)
# help(ksvm, help_type = "html")
m1 <- ksvm(x = X,
y = y,
scale = FALSE,
type = "C-svc",
kernel = "vanilladot",
C = 1)
m1
# Estrutura do objeto.
str(m1)
# Classe, arquitetura e métodos.
class(m1)
isS4(m1)
methods(class = class(m1))
# Vetores de suporte.
X[m1@SVindex, ]
TODO: extrair os parâmetros do hiperplano.
TODO: tabela de confusão.
TODO: classificação incorretas.
TODO: predição
```
# Classificação de cultivares de uva de vinho
Os dados são medidas de comprimento feitos em 100 folhas de 3 variedades
de uva: malbec, marlot e souvignon blanc. Os dados foram fornecidos pelo
pesquisador
[João Peterson Pereira Gardin](https://www.researchgate.net/profile/Joao_Gardin).
Os valores de área das folhas foram determinados por análise de imagem
das folhas digitalizadas por scanner usando o pacote
[EBImage](http://bioconductor.org/packages/release/bioc/html/EBImage.html).
```{r}
#-----------------------------------------------------------------------
# Dados hospedados na web.
url <- "http://www.leg.ufpr.br/~walmes/data/areafoliarUva.txt"
uva <- read.table(url, header = TRUE, sep = "\t",
stringsAsFactors = FALSE)
uva$cult <- factor(uva$cult)
uva$id <- NULL
# Comprimento da nervura lateral: média dos lados direito e esquerdo.
uva$nl <- with(uva, apply(cbind(nld, nle), 1, mean))
uva <- subset(uva, select = -c(nld, nle))
str(uva)
splom(uva[-(1:2)],
groups = uva$cult,
auto.key = TRUE,
cex = 0.2)
```
```{r, fig.height = 12}
splom(~uva[-(1:2)] | uva$cult,
cex = 0.2,
layout = c(1, NA))
```
## Pacote `e0171`
O código aperfeiçoado do material: <http://rischanlab.github.io/SVM.html>.
```{r}
library(e1071)
splom(iris[1:4],
groups = iris$Species)
xyplot(Petal.Length ~ Sepal.Length,
data = iris,
groups = Species,
auto.key = TRUE)
# Apenas duas espécies.
irisb <- droplevels(subset(iris,
Species != "virginica",
select = c(1, 3, 5)))
str(irisb)
xyplot(Petal.Length ~ Sepal.Length,
data = irisb,
groups = Species,
auto.key = TRUE)
# Exemplo com o Iris. Usando apenas duas variáveis.
x <- subset(irisb, select = -Species)
y <- irisb$Species
# Especificação com fórmula
svm0 <- svm(Species ~ Petal.Length + Sepal.Length,
data = irisb)
summary(svm0)
str(svm0)
# Quantidade e coordenadas dos pontos de suporte.
svm0$nSV
svm0$SV
# ATTENTION: O processamento dos dados é na escala padronizada!
plot(scale(Petal.Length) ~ scale(Sepal.Length),
data = irisb)
# points(x = svm0$SV[, "Sepal.Length"],
# y = svm0$SV[, "Petal.Length"],
# col = 2,
# pch = 19)
segments(x0 = 0,
y0 = 0,
x1 = svm0$SV[, "Sepal.Length"],
y1 = svm0$SV[, "Petal.Length"],
col = 2,
pch = 19)
abline(v = 0, h = 0, lty = 2, col = "gray")
# Usando matrizes e vetores.
svm1 <- svm(x, y)
summary(svm1)
# Predição dos valores.
pred <- predict(svm1, newdata = x)
# Performance do classificador.
table(pred, y)
# Gerando um grid no retângulo que contém os pontos.
grid <- with(irisb,
expand.grid(Sepal.Length = seq(min(Sepal.Length),
max(Sepal.Length),
l = 50),
Petal.Length = seq(min(Petal.Length),
max(Petal.Length),
l = 50)))
grid$y <- predict(svm1, newdata = grid)
# Exibindo a fronteira de classificação.
xyplot(Petal.Length ~ Sepal.Length,
data = grid,
groups = y,
pch = 4,
auto.key = TRUE) +
as.layer(xyplot(Petal.Length ~ Sepal.Length,
data = irisb,
pch = 19,
groups = Species))
#-----------------------------------------------------------------------
# Usando com todas as variáveis e 3 classes de espécie.
svm0 <- svm(Species ~ .,
data = iris)
summary(svm0)
# Predição.
yfit <- predict(svm0, newdata = iris)
# Tabela de confusão.
table(iris$Species, yfit)
#-----------------------------------------------------------------------
# Fazendo a tunagem.
x <- subset(iris, select = -Species)
y <- iris$Species
svm_tune <- tune(method = svm,
train.x = x,
train.y = y,
kernel = "radial",
ranges = list(cost = 10^(-1:2),
gamma = c(0.5, 1, 2)))
print(svm_tune)
# Usando os valores otimizados na validação cruzada.
svm1 <- svm(Species ~ .,
data = iris,
kernel = "radial",
cost = 1,
gamma = 0.5)
summary(svm1)
yfit <- predict(svm1, newdata = iris)
table(yfit, y)
```
## Pacote `kernlab`
```{r}
library(kernlab)
# help(ksvm, help_type = "html")
str(uva)
# Chamada com apenas duas classes. Simplificar para aprender.
da <- uva
levels(da$cult) <- c("malbec-merlot", "malbec-merlot", "sauvignonblanc")
table(da$cult)
m0 <- ksvm(cult ~ ., data = da)
m0
# Classe e funções disponíveis.
class(m0)
isS4(m0)
methods(class = class(m0))
# Número de vetores de suporte.
nSV(m0)
# Classficação nas observações de treino.
table(fitted(m0))
# Erro de classificação.
error(m0)
# Parâmetros.
param(m0)
# Só funciona para classificações binárias com duas preditoras.
# plot(m0)
splom(~da[, -1] | da$cult,
groups = fitted(m0),
auto.key = list(title = "Classificação"))
# Matriz de confusão.
ct <- table(fitted(m0), da$cult)
prop.table(ct)
# Percentual de acerto na classificação.
100 * sum(diag(ct))/sum(ct)
#-----------------------------------------------------------------------
# Parametrizando a chamada do método.
m0 <- ksvm(cult ~ .,
data = uva,
scaled = TRUE, # Padronizar com média 0 e variância 1.
type = "C-svc", # Emprego: classificação/regressão, etc.
C = 0.1, # Parâmetros do tipo de tarefa.
kernel = "rbfdot", # Função kernel.
kpar = list(sigma = 1), # Parâmetros da função kernel.
cross = 10) # Quantidade de folds para validação cruzada.
m0
# Matriz de confusão.
ct <- table(uva$cult, fitted(m0))
prop.table(ct)
# Gráfico de mosaico da matriz de confusão.
mosaicplot(ct,
color = brewer.pal(n = nlevels(uva$cult),
name = "Spectral"))
# Percentual de acerto na classificação.
100 * sum(diag(ct))/sum(ct)
#-----------------------------------------------------------------------
# Mais variações.
# Para diminuir o número de vetores de suporte.
m0 <- ksvm(cult ~ .,
data = uva,
scaled = TRUE,
type = "C-svc",
kernel = "rbfdot",
kpar = list(sigma = 0.01),
C = 1e7,
cross = 10)
m0
# Kernel linear (baunilha).
m0 <- ksvm(cult ~ .,
data = uva,
scaled = TRUE,
type = "nu-svc",
kernel = "vanilladot",
nu = 0.5,
cross = 10)
m0
# Percentual de acerto na classificação.
ct <- table(fitted(m0), uva$cult)
100 * sum(diag(ct))/sum(ct)
```
## Pacote `caret`
Incluir aplicações com o pacote `caret`, disponíveis nos tutoriais:
* http://dataaspirant.com/2017/01/19/support-vector-machine-classifier-implementation-r-caret-package/
* https://eight2late.wordpress.com/2017/02/07/a-gentle-introduction-to-support-vector-machines-using-r/
<!-- # Links úteis -->
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment