Skip to content

GitLab

machine-learning
machine-learning
Commit a897ae34 authored by Walmes Marques Zeviani's avatar Walmes Marques Zeviani
Browse files
Options

Melhora o tutorial de fundamentos de SVM.

parent e7ba2d59
Hide whitespace changes
Inline Side-by-side
Showing with 93 additions and 464 deletions
tutorials/05-support-vector-machine.Rmd
View file @ a897ae34
......@@ -22,6 +22,13 @@ library(lattice)
library(latticeExtra)
```
Nessa primeira sessão será feita uma exposição dos componentes básicos
relacionados ao algorítmo de máquinas de vetores de suposte para
classificação. Na primeira parte será ilistrado o caso linearmente
separável e depois o caso não linearmente separável. Será usando o
pacote `quadprog` para resolver o problema de otimização quadrática com
inequações lineares.
## Simulando um caso linearmente separável
O código abaixo apenas ilustra um caso simples de classificação em duas
......@@ -38,11 +45,11 @@ 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"))
da$y <- with(da, ifelse(0.1 + x1 - x2 < 0, "cyan", "orange"))
# Grafico com as categorias com cores diferentes.
with(da, plot(x2 ~ x1, col = y, asp = 1))
abline(a = 0.1, b = 1, col = "magenta")
with(da, plot(x2 ~ x1, col = y, asp = 1, pch = 19))
abline(a = 0.1, b = 1, col = "purple", lwd = 2)
```
Uma solução construída usando o MATLAB:
......@@ -50,21 +57,23 @@ Uma solução construída usando o MATLAB:
## Transformando em não linearmente separável
Para gerar uma sepação não linear, um ruído branco será somado aos
valores observados das variáveis métricas de cada registro.
Para gerar uma sepação não linear, um ruído branco (perturbação) será
somado aos valores observados das variáveis métricas de cada
registro. Dessa forma, os dados serão não linearmente separaveis, mais
condizente com a maioria dos problemas de classificação reais.
```{r}
#-----------------------------------------------------------------------
# Versão com um jitter nos pontos para dar uma misturada.
# Versão com uma perturbação nos pontos para dar uma bagunçada.
set.seed(1234)
a <- 0.2
a <- 0.2 # Semi-amplitude da perturbação.
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")
with(db, plot(x2 ~ x1, col = y, asp = 1, pch = 19))
abline(a = 0.1, b = 1, col = "purple")
```
## Otimização quadrática com restrições
......@@ -84,28 +93,34 @@ foi apresentada no stackoverflow:
<https://stackoverflow.com/questions/33088363/how-to-use-r-package-quadprog-to-solve-svm>.
Aqui o código foi revisado e ampliado.
Uma abordagem usando o Julia foi feita pelo Professor Abel Siqueira
(Departamento de Matemática - UFPR):
<http://abelsiqueira.github.io/disciplinas/cm116/2018/SVM.html>.
```{r}
# Pacote para problemas de otimização convexa com restrições lineares.
library(quadprog)
ls("package:quadprog")
# Acesso à documentação.
# help(package = "quadprog", h = "html")
# help(package = "quadprog", help_type = "html")
# Aplicando para os dados de duas especies de iris com as variáveis de
# 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.
# Matriz de variáveis preditoras normalizadas.
X <- scale(as.matrix(da[, 1:2]))
# Para deslocar os dados da coordenada (0, 0).
X[, 1] <- X[, 1] + 1
X[, 2] <- X[, 2] - 1
# Cabeça e cauda das matrizes.
head(cbind(X, y))
tail(cbind(X, y))
......@@ -115,7 +130,7 @@ plot(X, col = y + 3)
grid()
#--------------------------------------------
# Desenvolvendo a solução com o quadprog.
# Desenvolvendo a solução com o `quadprog`.
# Número de observações.
N <- nrow(X)
......@@ -125,7 +140,8 @@ N
n_d <- ncol(X)
n_d
# Valor do parâmetro de regularização.
# Valor do hiperparâmetro de regularização (parâmetro que penaliza a
# soma das distâncias da margem das observações classificadas errado).
C <- 1
```
......@@ -151,9 +167,12 @@ d <- c(0,
rep(0, n_d),
rep(-C, N))
# Número de variáveis na otimização.
length(d)
# Valor positivo pequeno mas não nulo para fazer a matriz D ser positiva
# definida (exigência da quadprog()).
eps <- 1e-10
eps <- 1e-14
D <- diag(c(eps,
rep(1, n_d),
rep(eps, N)))
......@@ -192,12 +211,18 @@ b_0 <- c(rep(0, N),
length(b_0)
# Passando os componentes para o otimizador.
results <- solve.QP(D, d, A, b_0)
results <- solve.QP(Dmat = D,
dvec = d,
Amat = A,
bvec = b_0)
# Extrai os parâmetros otimizados, i.e., [beta_0, beta, zeta].
b_optim <- results$solution
b_optim
b_optim[1:3] # [beta_0, beta].
b_optim[-(1:3)] # [zeta], se zeta_i != 0 então é ponto de suporte.
# Hiperplano de separação.
beta_0 <- b_optim[1]
beta_0
......@@ -205,9 +230,11 @@ beta_0
beta <- b_optim[1 + (1:n_d)]
beta
# 1/sum(beta^2)
# Tamanho da margem.
m <- 1/sum(beta^2)
m
# Coeficientes.
# Coeficientes associados a cada ponto.
zeta <- b_optim[(n_d + 1) + (1:N)]
zeta
......@@ -216,9 +243,9 @@ C * sum(zeta) + 0.5 * sum(beta^2)
results$value
# Vetores de suporte.
sp <- rbind(X[which(zeta > 0), ])
sp <- rbind(X[which(zeta != 0), ])
sp
nrow(sp)
c(nrow(sp), nrow(unique(sp)))
# Predição a partir do modelo ajustado.
y_pred <- sign(apply(X,
......@@ -235,8 +262,8 @@ 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),
seq(min(X[, 1]), max(X[, 1]), length.out = 61),
seq(min(X[, 2]), max(X[, 2]), length.out = 61),
KEEP.OUT.ATTRS = FALSE))
names(grid) <- names(da)[1:2]
grid$y_pred <- sign(apply(as.matrix(grid),
......@@ -252,12 +279,22 @@ plot(X,
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")
col = c("orange", "cyan")[0.5 * grid$y_pred + 1.5],
pch = 3,
cex = 0.5)
for (d in c(-1, 0, 1)) {
abline(a = -beta_0/beta[2] + beta[2] * m * d,
b = -beta[-2]/beta[2],
lty = 2,
col = "black")
}
abline(v = 0, h = 0, lty = 2, col = "gray")
arrows(x0 = 0, x1 = c(-1, 1) * beta[1] * m,
y0 = 0, y1 = c(-1, 1) * beta[2] * m,
lty = 1,
lwd = 2,
length = 0.05,
col = "green4")
arrows(x0 = 0,
y0 = 0,
x1 = sp[, 1],
......@@ -269,7 +306,7 @@ arrows(x0 = 0,
if (length(i)) {
points(rbind(X[i, ]),
pch = 19,
cex = 0.8,
cex = 0.6,
col = "yellow")
}
```
......@@ -294,9 +331,10 @@ C++.
```{r}
library(e1071)
ls("package:e1071")
# Acessa a documentação da função.
# help(svm, help_type = "html")
# help(svm, package = "e1071", help_type = "html")
m0 <- svm(x = X,
y = y,
......@@ -310,6 +348,10 @@ m0
class(m0)
methods(class = class(m0))
# Abre a função predict.
# getS3method(f = "predict", class = "svm")
str(m0$decision.values)
# Tabela de confusão.
table(y, predict(m0))
......@@ -331,12 +373,9 @@ plot(X,
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")
col = c("orange", "cyan")[as.integer(grid$y_svm)],
pch = 3,
cex = 0.5)
arrows(x0 = 0,
y0 = 0,
x1 = sp[, 1],
......@@ -348,7 +387,7 @@ arrows(x0 = 0,
if (length(i)) {
points(rbind(X[i, ]),
pch = 19,
cex = 0.8,
cex = 0.6,
col = "yellow")
}
```
......@@ -359,15 +398,22 @@ O pacote [`kernlab`](https://cran.r-project.org/web/packages/kernlab)
possui a função `ksvm()` que também se comunica com a `libsvm`. No
entanto, a arquitetura do pacote `kernlab` é S4.
A diferença é que o `kernlab` permite empregar o truque do kernel porque
possui várias funções kernel implementadas. Funções kernel serão
consideradas adiante.
```{r}
library(kernlab)
# Funções kernel.
grep(ls("package:kernlab"), pattern = "dot$", value = TRUE)
# help(ksvm, help_type = "html")
m1 <- ksvm(x = X,
y = y,
scale = FALSE,
type = "C-svc",
kernel = "vanilladot",
kernel = "vanilladot", # Sem função kernel.
C = 1)
m1
......@@ -403,8 +449,9 @@ plot(X,
cex = 1.2,
pch = 19)
points(grid[, 1:2],
col = (grid$y_ksvm + 3),
pch = 3)
col = c("orange", "cyan")[as.integer(grid$y_svm)],
pch = 3,
cex = 0.5)
arrows(x0 = 0,
y0 = 0,
x1 = sp[, 1],
......@@ -416,425 +463,15 @@ arrows(x0 = 0,
if (length(i)) {
points(rbind(X[i, ]),
pch = 19,
cex = 0.8,
cex = 0.6,
col = "yellow")
}
```
<!-- TODO: utilizar outras funções kernel. -->
# 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))
```
```{r, eval = FALSE, include = FALSE}
# Testando componentes principais.
u <- scale(uva[, -1])
pr <- princomp(x = u)
screeplot(pr, type = "lines")
biplot(pr)
plot(pr$scores[, 1:2],
col = as.integer(uva$cult))
abline(v = 0, h = 0, lty = 2)
```
## 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 = 1e5,
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 `e1071`
```{r}
library(e1071)
# Especificação com fórmula.
m1 <- svm(cult ~ .,
data = uva)
summary(m1)
str(m1)
# Quantidade e coordenadas dos pontos de suporte.
m1$tot.nSV
m1$nSV
head(m1$SV)
#-----------------------------------------------------------------------
# Fazendo a tunagem com a grid search.
x <- as.matrix(subset(uva, select = -cult))
y <- uva$cult
tune <- tune(method = svm,
train.x = x,
train.y = y,
kernel = "radial",
ranges = list(cost = 10^seq(-2, 5, l = 10),
gamma = 2^seq(-3, 3, l = 10)))
print(tune)
str(tune)
# A superfície do erro.
levelplot(error ~ log10(cost) + log2(gamma),
data = tune$performances,
contour = TRUE) +
layer(panel.abline(v = log10(cost), h = log2(gamma), lty = 2),
data = tune$best.model)
# Usando os valores otimizados na validação cruzada.
m2 <- svm(cult ~ .,
data = uva,
kernel = "radial",
cost = tune$best.model$cost,
gamma = tune$best.model$gamma)
summary(m2)
yfit <- predict(m2)
ct <- table(yfit, y)
sum(diag(ct))/sum(ct)
```
O código abaixo foi aperfeiçoado do material
<http://rischanlab.github.io/SVM.html> e trabalha os dados de especies
de iris em `iris`.
```{r}
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)
#
# # 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 = 101),
# Petal.Length = seq(min(Petal.Length),
# max(Petal.Length),
# l = 101)))
# grid$y <- predict(svm1, newdata = grid)
#
# # Exibindo a fronteira de classificação.
# xyplot(Petal.Length ~ Sepal.Length,
# data = grid,
# groups = y,
# pch = 3,
# aspect = 1,
# 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 e avaliando
# diferentes funções kernel.
# Nomes das funções kernel.
ker <- c("linear", "polynomial", "radial", "sigmoid")
svm0 <- sapply(ker,
simplify = FALSE,
FUN = function(k) {
svm(Species ~ Petal.Length + Sepal.Length,
data = iris,
kernel = k)
})
lapply(svm0, summary)
# Tabelas de confusão.
lapply(svm0,
FUN = function(model) {
ct <- table(iris$Species, predict(model))
cat("Acertos:", 100 * sum(diag(ct))/sum(ct), "\n")
return(ct)
})
# Gerando um grid no retângulo que contém os pontos.
grid <- with(iris,
expand.grid(Sepal.Length = seq(min(Sepal.Length),
max(Sepal.Length),
l = 51),
Petal.Length = seq(min(Petal.Length),
max(Petal.Length),
l = 51)))
y <- lapply(svm0, FUN = predict, newdata = grid)
y <- lapply(y, as.data.frame)
names(y) <- ker
y <- plyr::ldply(y)
names(y) <- c("kernel", "y")
grid <- cbind(grid, y)
str(grid)
# Exibindo a fronteira de classificação.
xyplot(Petal.Length ~ Sepal.Length | kernel,
data = grid,
groups = y,
pch = 3,
as.table = TRUE,
aspect = 1,
auto.key = TRUE) +
as.layer(xyplot(Petal.Length ~ Sepal.Length,
data = iris,
pch = 19,
groups = Species))
```
## Pacote `caret`
Adaptação feita baseada no material <http://dataaspirant.com/2017/01/19/support-vector-machine-classifier-implementation-r-caret-package/>.
```{r}
library(caret)
set.seed(1234)
# Especificação da validação cruzada.
trctrl <- trainControl(method = "repeatedcv",
number = 10,
repeats = 3)
# Ajuste.
svm_Linear <- train(cult ~ .,
data = uva,