Atualiza conteúdo de SVM.

parent 5f316d80
......@@ -24,6 +24,12 @@ library(latticeExtra)
## Simulando um caso linearmente separável
O código abaixo apenas ilustra um caso simples de classificação em duas
categorias a partir de duas variáveis métricas. Os dados simulados tem
uma estrutura bastante regular para facilitar o desenvolvimento e
verificação da implementação de classificadores baseados em vetores de
suporte.
```{r}
#-----------------------------------------------------------------------
# Gerando um conjunto de dados artifical.
......@@ -34,6 +40,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"))
# Grafico com as categorias com cores diferentes.
with(da, plot(x2 ~ x1, col = y, asp = 1))
abline(a = 0.1, b = 1, col = "magenta")
```
......@@ -43,12 +50,15 @@ 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.
```{r}
#-----------------------------------------------------------------------
# Versão com um jitter nos pontos para dar uma misturada.
set.seed(1234)
a <- 0.1
a <- 0.2
db <- transform(da,
x1 = x1 + runif(nrow(da), -a, a),
x2 = x2 + runif(nrow(da), -a, a))
......@@ -59,6 +69,10 @@ abline(a = 0.1, b = 1, col = "magenta")
## Otimização quadrática com restrições
Nessa seção do tutorial, o problema de otimização convexa com restrições
na forma de inequações lineares será resolvido utilizando o pacote
`quadprog`.
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
......@@ -68,6 +82,7 @@ otimização convexa com o pacote
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>.
Aqui o código foi revisado e ampliado.
```{r}
# Pacote para problemas de otimização convexa com restrições lineares.
......@@ -97,19 +112,39 @@ tail(cbind(X, y))
# Diagrama de dispersão com indicação das classes.
plot(X, col = y + 3)
grid()
#--------------------------------------------
# Desenvolvendo a solução com o quadprog.
# Número de observações.
N <- nrow(X)
N
# Número de dimensões.
n_d <- ncol(X)
n_d
# Valor do parâmetro de regularização.
C <- 1
```
A quantidade a ser minimizada é
$$
\min_{\beta_0 \in \mathbb{R},
\beta \in \mathbb{R}^p,
\xi \in \mathbb{R}^n}\,
C \sum_{i=1}^{n} \xi_i + \frac{1}{2} \beta^\top \beta,
$$
considerando as restrições
$$
\begin{align*}
\xi_i \geq 0\,\, \forall i\\
\xi_i + \beta_0 y_i + \beta^\top x_i y_i \geq 1\,\, \forall i.
\end{align*}
$$
```{r}
# Vetor que multiplica o vetor de parâmetros
# b = [\beta_0, \beta, \zeta].
d <- c(0,
......@@ -159,9 +194,6 @@ 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
......@@ -179,6 +211,10 @@ beta
zeta <- b_optim[(n_d + 1) + (1:N)]
zeta
# Valor da função objetivo na convergência.
C * sum(zeta) + 0.5 * sum(beta^2)
results$value
# Vetores de suporte.
sp <- rbind(X[which(zeta > 0), ])
sp
......@@ -238,10 +274,6 @@ if (length(i)) {
}
```
# 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
......@@ -323,6 +355,10 @@ if (length(i)) {
## Pacote `kernlab`
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.
```{r}
library(kernlab)
......@@ -343,15 +379,50 @@ class(m1)
isS4(m1)
methods(class = class(m1))
# Valor da função objetivo.
m1@obj
# Vetores de suporte.
X[m1@SVindex, ]
# X[m1@SVindex, ]
sp <- m1@xmatrix[[1]]
sp
# Tabela de confusão.
table(m1@ymatrix, m1@fitted)
# Classificações incorretas.
i <- which(m1@ymatrix * m1@fitted < 0)
# Fazendo a predição.
grid$y_ksvm <- predict(m1, newdata = grid[, 1:2])
TODO: extrair os parâmetros do hiperplano.
TODO: tabela de confusão.
TODO: classificação incorretas.
TODO: predição
# 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_ksvm + 3),
pch = 3)
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")
}
```
<!-- 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
......@@ -388,130 +459,18 @@ splom(~uva[-(1:2)] | uva$cult,
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)
```{r, eval = FALSE, include = FALSE}
# Testando componentes principais.
u <- scale(uva[, -1])
pr <- princomp(x = u)
xyplot(Petal.Length ~ Sepal.Length,
data = irisb,
groups = Species,
auto.key = TRUE)
screeplot(pr, type = "lines")
# Exemplo com o Iris. Usando apenas duas variáveis.
x <- subset(irisb, select = -Species)
y <- irisb$Species
biplot(pr)
# 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)
plot(pr$scores[, 1:2],
col = as.integer(uva$cult))
abline(v = 0, h = 0, lty = 2)
```
## Pacote `kernlab`
......@@ -596,7 +555,7 @@ m0 <- ksvm(cult ~ .,
type = "C-svc",
kernel = "rbfdot",
kpar = list(sigma = 0.01),
C = 1e7,
C = 1e5,
cross = 10)
m0
......@@ -615,12 +574,337 @@ 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/>.
Incluir aplicações com o pacote `caret`, disponíveis nos tutoriais:
```{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,
method = "svmLinear",
trControl = trctrl,
preProcess = c("center", "scale"),
tuneLength = 10)
# Resultado do ajuste.
svm_Linear
str(svm_Linear)
# Mostra que foi feito a chamada da kernlab::ksvm().
svm_Linear$finalModel
# Matriz de confusão.
confusionMatrix(uva$cult, predict(svm_Linear))
#-----------------------------------------------------------------------
# Tunando.
grid <- expand.grid(C = c(0, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75,
1, 1.25, 1.5, 1.75, 2, 5))
svm_Linear_Grid <- train(cult ~ .,
data = uva,
method = "svmLinear",
trControl = trctrl,
preProcess = c("center", "scale"),
tuneGrid = grid,
tuneLength = 10)
* 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/
svm_Linear_Grid
plot(svm_Linear_Grid)
#-----------------------------------------------------------------------
# Usando kernel radial.
<!-- # Links úteis -->
svm_Radial <- train(Species ~ .,
data = iris,
method = "svmRadial",
trControl = trctrl,
preProcess = c("center", "scale"),
tuneLength = 10)
svm_Radial
plot(svm_Radial)
```
# Funções kernel
```{r}
# help(dots, package = "kernlab", h = "html")
# Revelando o código de uma função kernel.
rbfkernel <- rbfdot(sigma = 0.1)
rbfkernel
# Parâmetros da função.
kpar(rbfkernel)
# Mostra o corpo da função.
body(rbfkernel)
# Função kernel definida pelo usuário (tirado da documentação).
k <- function(x, y) {
(sum(x * y) + 1) * exp(-0.001 * sum((x - y)^2))
}
class(k) <- "kernel"
# Avaliando a função kernel.
x <- runif(10)
y <- runif(10)
k(x, y)
# Usando a função kernel definida pelo usuário.
m0 <- ksvm(Species ~ .,
data = iris,
kernel = k,
C = 1,
cross = 10)
m0
#-----------------------------------------------------------------------
# Expansão de características por meio de funções kernel.
# help(kernelMatrix, h = "html")
# Duas características.
x1 <- cbind(seq(0, 10, by = 1))
x2 <- cbind(seq(-1, 1, length.out = 11))
cbind(x1, x2)
# Kernel linear.
body(vanilladot())