1 Classificador de vetores de suporte
2 Implementações de máquinas de vetores de suporte
- 2.1 Pacote e1071
- 2.2 Pacote kernlab
3 Funções kernel
- 3.1 Explorando as funções kernel
- 3.2 Ajustes com diferentes funções kernel
  - 3.2.1 Linear

1 Classificador de vetores de suporte

rm(list = objects())
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.

1.1 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.

#-----------------------------------------------------------------------
# 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, "cyan", "orange"))

# Grafico com as categorias com cores diferentes.
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: https://oceandatamining.sciencesconf.org/conference/oceandatamining/TP_Linear_SVM.pdf.

1.2 Transformando em não linearmente separável

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.

#-----------------------------------------------------------------------
# Versão com uma perturbação nos pontos para dar uma bagunçada.

set.seed(1234)
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, pch = 19))
abline(a = 0.1, b = 1, col = "purple")

1.3 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 otimização convexa com o pacote 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. 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.

# Pacote para problemas de otimização convexa com restrições lineares.
library(quadprog)
ls("package:quadprog")

## [1] "solve.QP"         "solve.QP.compact"

# Acesso à documentação.
# help(package = "quadprog", help_type = "html")

# Matriz de diagramas de dispersão.
splom(iris[, 1:4], groups = iris$Species)

# Aplicando para os dados de duas especies de `iris` com as variáveis de
# comprimento apenas.
da <- droplevels(subset(iris,
                        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 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))

##    Sepal.Length Petal.Length  y
## 51    2.1134002    -1.249522 -1
## 52    1.2081968    -1.491776 -1
## 53    1.9625330    -1.007268 -1
## 54   -0.1496083    -2.097412 -1
## 55    1.3590640    -1.370649 -1
## 56    0.1521261    -1.491776 -1

tail(cbind(X, y))

##     Sepal.Length Petal.Length y
## 145    1.6607985   -0.0382501 1
## 146    1.6607985   -0.6438861 1
## 147    1.0573295   -0.8861404 1
## 148    1.3590640   -0.6438861 1
## 149    0.9064623   -0.4016317 1
## 150    0.4538606   -0.7650132 1

# 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

## [1] 100

# Número de dimensões.
n_d <- ncol(X)
n_d

## [1] 2

# 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

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*} \]

# Vetor que multiplica o vetor de parâmetros
# b = [\beta_0, \beta, \zeta].
d <- c(0,
       rep(0, n_d),
       rep(-C, N))

# Número de variáveis na otimização.
length(d)

## [1] 103

# Valor positivo pequeno mas não nulo para fazer a matriz D ser positiva
# definida (exigência da quadprog()).
eps <- 1e-14
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)

## [1] 100 103

# 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)

## [1] 100 103

# 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)

## [1] 103 200

# image(A, asp = 1)

b_0 <- c(rep(0, N),
         rep(1, N))
length(b_0)

## [1] 200

# Passando os componentes para o otimizador.
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

##   [1]  3.847948e+00 -6.034689e-01  3.149307e+00  0.000000e+00  0.000000e+00
##   [6]  4.914251e-01  0.000000e+00  0.000000e+00  5.808260e-02  2.747539e-01
##  [11]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [16]  0.000000e+00  4.568412e-01  0.000000e+00  0.000000e+00  1.491263e-01
##  [21]  0.000000e+00  0.000000e+00  0.000000e+00  1.020395e+00  0.000000e+00
##  [26]  1.037687e+00  4.568412e-01  0.000000e+00  0.000000e+00  2.010021e-01
##  [31]  1.054979e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [36]  0.000000e+00  2.073752e+00  3.312137e-01  0.000000e+00  0.000000e+00
##  [41]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  7.537453e-02
##  [46]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [51]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [56]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.213568e+00
##  [61]  0.000000e+00  0.000000e+00  0.000000e+00  3.814667e-01  0.000000e+00
##  [66]  0.000000e+00  3.458388e-02  0.000000e+00  0.000000e+00  0.000000e+00
##  [71]  0.000000e+00  0.000000e+00  3.077150e-01  0.000000e+00  3.250069e-01
##  [76]  0.000000e+00  9.623127e-01  0.000000e+00  0.000000e+00  1.252736e+00
##  [81]  7.802253e-01  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [86]  0.000000e+00  1.993793e-01  6.952686e-16  0.000000e+00  0.000000e+00
##  [91]  0.000000e+00  1.070648e+00  0.000000e+00  0.000000e+00  7.456415e-01
##  [96]  0.000000e+00  0.000000e+00  0.000000e+00  1.820874e-01  5.808460e-01
## [101]  0.000000e+00  0.000000e+00  0.000000e+00

b_optim[1:3]    # [beta_0, beta].

## [1]  3.8479478 -0.6034689  3.1493070

b_optim[-(1:3)] # [zeta], se zeta_i != 0 então é ponto de suporte.

##   [1] 0.000000e+00 0.000000e+00 4.914251e-01 0.000000e+00 0.000000e+00
##   [6] 5.808260e-02 2.747539e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [11] 0.000000e+00 0.000000e+00 0.000000e+00 4.568412e-01 0.000000e+00
##  [16] 0.000000e+00 1.491263e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [21] 1.020395e+00 0.000000e+00 1.037687e+00 4.568412e-01 0.000000e+00
##  [26] 0.000000e+00 2.010021e-01 1.054979e+00 0.000000e+00 0.000000e+00
##  [31] 0.000000e+00 0.000000e+00 0.000000e+00 2.073752e+00 3.312137e-01
##  [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [41] 0.000000e+00 7.537453e-02 0.000000e+00 0.000000e+00 0.000000e+00
##  [46] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [51] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [56] 0.000000e+00 1.213568e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [61] 3.814667e-01 0.000000e+00 0.000000e+00 3.458388e-02 0.000000e+00
##  [66] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 3.077150e-01
##  [71] 0.000000e+00 3.250069e-01 0.000000e+00 9.623127e-01 0.000000e+00
##  [76] 0.000000e+00 1.252736e+00 7.802253e-01 0.000000e+00 0.000000e+00
##  [81] 0.000000e+00 0.000000e+00 0.000000e+00 1.993793e-01 6.952686e-16
##  [86] 0.000000e+00 0.000000e+00 0.000000e+00 1.070648e+00 0.000000e+00
##  [91] 0.000000e+00 7.456415e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [96] 1.820874e-01 5.808460e-01 0.000000e+00 0.000000e+00 0.000000e+00

# Hiperplano de separação.
beta_0 <- b_optim[1]
beta_0

## [1] 3.847948

beta <- b_optim[1 + (1:n_d)]
beta

## [1] -0.6034689  3.1493070

# Tamanho da margem.
m <- 1/sum(beta^2)
m

## [1] 0.09725442

# Coeficientes associados a cada ponto.
zeta <- b_optim[(n_d + 1) + (1:N)]
zeta

##   [1] 0.000000e+00 0.000000e+00 4.914251e-01 0.000000e+00 0.000000e+00
##   [6] 5.808260e-02 2.747539e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [11] 0.000000e+00 0.000000e+00 0.000000e+00 4.568412e-01 0.000000e+00
##  [16] 0.000000e+00 1.491263e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [21] 1.020395e+00 0.000000e+00 1.037687e+00 4.568412e-01 0.000000e+00
##  [26] 0.000000e+00 2.010021e-01 1.054979e+00 0.000000e+00 0.000000e+00
##  [31] 0.000000e+00 0.000000e+00 0.000000e+00 2.073752e+00 3.312137e-01
##  [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [41] 0.000000e+00 7.537453e-02 0.000000e+00 0.000000e+00 0.000000e+00
##  [46] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [51] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [56] 0.000000e+00 1.213568e+00 0.000000e+00 0.000000e+00 0.000000e+00
##  [61] 3.814667e-01 0.000000e+00 0.000000e+00 3.458388e-02 0.000000e+00
##  [66] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 3.077150e-01
##  [71] 0.000000e+00 3.250069e-01 0.000000e+00 9.623127e-01 0.000000e+00
##  [76] 0.000000e+00 1.252736e+00 7.802253e-01 0.000000e+00 0.000000e+00
##  [81] 0.000000e+00 0.000000e+00 0.000000e+00 1.993793e-01 6.952686e-16
##  [86] 0.000000e+00 0.000000e+00 0.000000e+00 1.070648e+00 0.000000e+00
##  [91] 0.000000e+00 7.456415e-01 0.000000e+00 0.000000e+00 0.000000e+00
##  [96] 1.820874e-01 5.808460e-01 0.000000e+00 0.000000e+00 0.000000e+00

# Valor da função objetivo na convergência.
C * sum(zeta) + 0.5 * sum(beta^2)

## [1] 20.85885

results$value

## [1] 20.85885

# Vetores de suporte.
sp <- rbind(X[which(zeta != 0), ])
sp

##     Sepal.Length Petal.Length
## 53   1.962532967   -1.0072676
## 56   0.152126133   -1.4917764
## 57   1.057329550   -1.2495220
## 64   0.755595077   -1.2495220
## 67   0.001258897   -1.4917764
## 71   0.453860605   -1.1283948
## 73   1.057329550   -1.0072676
## 74   0.755595077   -1.2495220
## 77   1.811665730   -1.1283948
## 78   1.660798494   -0.8861404
## 84   0.604727841   -0.7650132
## 85  -0.300475576   -1.4917764
## 92   0.755595077   -1.3706492
## 107 -1.054811756   -1.4917764
## 111  1.359064022   -0.7650132
## 114  0.152126133   -0.8861404
## 120  0.604727841   -0.8861404
## 122  0.001258897   -1.0072676
## 124  1.057329550   -1.0072676
## 127  0.906462314   -1.1283948
## 128  0.755595077   -1.0072676
## 134  1.057329550   -0.7650132
## 135  0.755595077   -0.1593773
## 139  0.604727841   -1.1283948
## 142  1.962532967   -0.7650132
## 146  1.660798494   -0.6438861
## 147  1.057329550   -0.8861404

c(nrow(sp), nrow(unique(sp)))

## [1] 27 25

# 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)

##     y_pred
## y    -1  1
##   -1 46  4
##   1   3 47

# Classificações incorretas.
i <- which(y * y_pred < 0)
length(i)

## [1] 7

# Fazendo um grid para predição.
grid <- with(da,
             expand.grid(
                 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),
                          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 = 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],
       y1 = sp[, 2],
       pch = 4,
       cex = 2,
       col = "magenta",
       length = 0.05)
if (length(i)) {
    points(rbind(X[i, ]),
           pch = 19,
           cex = 0.6,
           col = "yellow")
}

2 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.

2.1 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 possui a função svm() que dá acesso a implementação disponível na libsvm desenvolvida em C++. A e1071::svm() tem 4 opções para funções kernel que serão comentadas posteriormente. Consulte a documentação para mais detalhes.

library(e1071)
ls("package:e1071")

##  [1] "allShortestPaths"      "bclust"               
##  [3] "best.nnet"             "best.randomForest"    
##  [5] "best.rpart"            "best.svm"             
##  [7] "best.tune"             "bincombinations"      
##  [9] "bootstrap.lca"         "centers.bclust"       
## [11] "classAgreement"        "clusters.bclust"      
## [13] "cmeans"                "compareMatchedClasses"
## [15] "countpattern"          "cshell"               
## [17] "d2sigmoid"             "ddiscrete"            
## [19] "dsigmoid"              "element"              
## [21] "extractPath"           "fclustIndex"          
## [23] "hamming.distance"      "hamming.window"       
## [25] "hanning.window"        "hclust.bclust"        
## [27] "hsv_palette"           "ica"                  
## [29] "impute"                "interpolate"          
## [31] "kurtosis"              "lca"                  
## [33] "matchClasses"          "matchControls"        
## [35] "moment"                "naiveBayes"           
## [37] "pdiscrete"             "permutations"         
## [39] "probplot"              "qdiscrete"            
## [41] "rbridge"               "rdiscrete"            
## [43] "read.matrix.csr"       "rectangle.window"     
## [45] "rwiener"               "sigmoid"              
## [47] "skewness"              "stft"                 
## [49] "svm"                   "tune"                 
## [51] "tune.control"          "tune.knn"             
## [53] "tune.nnet"             "tune.randomForest"    
## [55] "tune.rpart"            "tune.svm"             
## [57] "write.matrix.csr"      "write.svm"

# Acessa a documentação da função.
# help(svm, package = "e1071", help_type = "html")

# Aplica SVM para classificação.
m0 <- svm(x = X,
          y = y,
          scale = FALSE,
          kernel = "linear",
          type = "C-classification",
          cost = 1)
m0

## 
## Call:
## svm.default(x = X, y = y, scale = FALSE, type = "C-classification", 
##     kernel = "linear", cost = 1)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.5 
## 
## Number of Support Vectors:  26

# Classe e métodos.
class(m0)

## [1] "svm"

methods(class = class(m0))

## [1] plot    predict print   summary
## see '?methods' for accessing help and source code

# Abre a função predict.
# getS3method(f = "predict", class = "svm")
str(m0$decision.values)

##  num [1:100, 1] 1.328 1.545 0.474 2.633 1.254 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:100] "51" "52" "53" "54" ...
##   ..$ : chr "-1/1"

# Tabela de confusão.
table(y, predict(m0))

##     
## y    -1  1
##   -1 46  4
##   1   3 47

# Classificações incorretas.
i <- which(y * as.integer(as.character(predict(m0))) < 0)
length(i)

## [1] 7

# Vetores de suporte.
sp <- m0$SV
sp

##     Sepal.Length Petal.Length
## 53   1.962532967   -1.0072676
## 56   0.152126133   -1.4917764
## 57   1.057329550   -1.2495220
## 64   0.755595077   -1.2495220
## 67   0.001258897   -1.4917764
## 71   0.453860605   -1.1283948
## 73   1.057329550   -1.0072676
## 74   0.755595077   -1.2495220
## 77   1.811665730   -1.1283948
## 78   1.660798494   -0.8861404
## 84   0.604727841   -0.7650132
## 85  -0.300475576   -1.4917764
## 92   0.755595077   -1.3706492
## 107 -1.054811756   -1.4917764
## 111  1.359064022   -0.7650132
## 114  0.152126133   -0.8861404
## 120  0.604727841   -0.8861404
## 122  0.001258897   -1.0072676
## 124  1.057329550   -1.0072676
## 127  0.906462314   -1.1283948
## 128  0.755595077   -1.0072676
## 134  1.057329550   -0.7650132
## 139  0.604727841   -1.1283948
## 142  1.962532967   -0.7650132
## 146  1.660798494   -0.6438861
## 147  1.057329550   -0.8861404

# 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 = c("orange", "cyan")[as.integer(grid$y_svm)],
       pch = 3,
       cex = 0.5)
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.6,
           col = "yellow")
}

2.2 Pacote `kernlab`

O pacote 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 explorar mais o truque do kernel porque possui várias funções kernel implementadas além de permitir que o usuário defina suas próprias funções kernel. Funções kernel serão consideradas adiante.

library(kernlab)

# Acesso à documentação.
# help(ksvm, package = "kernlab", help_type = "html")

# Funções kernel.
grep(ls("package:kernlab"), pattern = "dot$", value = TRUE)

## [1] "anovadot"   "besseldot"  "laplacedot" "polydot"    "rbfdot"    
## [6] "splinedot"  "stringdot"  "tanhdot"    "vanilladot"

# help(ksvm, help_type = "html")
m1 <- ksvm(x = X,
           y = y,
           scale = FALSE,
           type = "C-svc",
           kernel = "vanilladot", # Kernel linear.
           C = 1)

##  Setting default kernel parameters

m1

## Support Vector Machine object of class "ksvm" 
## 
## SV type: C-svc  (classification) 
##  parameter : cost C = 1 
## 
## Linear (vanilla) kernel function. 
## 
## Number of Support Vectors : 26 
## 
## Objective Function Value : -20.8588 
## Training error : 0.07

# Estrutura do objeto.
# str(m1)

# Classe, arquitetura e métodos.
class(m1)

## [1] "ksvm"
## attr(,"package")
## [1] "kernlab"

isS4(m1)

## [1] TRUE

methods(class = class(m1))

##  [1] alphaindex alpha      b          coef       cross      error     
##  [7] fitted     kcall      kernelf    kpar       lev        nSV       
## [13] obj        param      plot       predict    prior      prob.model
## [19] scaling    show       SVindex    type       xmatrix    ymatrix   
## see '?methods' for accessing help and source code

# Valor da função objetivo.
m1@obj

## [1] -20.85885

# Vetores de suporte.
# X[m1@SVindex, ]
sp <- m1@xmatrix[[1]]
sp

##     Sepal.Length Petal.Length
## 53   1.962532967   -1.0072676
## 56   0.152126133   -1.4917764
## 57   1.057329550   -1.2495220
## 64   0.755595077   -1.2495220
## 67   0.001258897   -1.4917764
## 71   0.453860605   -1.1283948
## 73   1.057329550   -1.0072676
## 74   0.755595077   -1.2495220
## 77   1.811665730   -1.1283948
## 78   1.660798494   -0.8861404
## 84   0.604727841   -0.7650132
## 85  -0.300475576   -1.4917764
## 92   0.755595077   -1.3706492
## 107 -1.054811756   -1.4917764
## 111  1.359064022   -0.7650132
## 114  0.152126133   -0.8861404
## 120  0.604727841   -0.8861404
## 122  0.001258897   -1.0072676
## 124  1.057329550   -1.0072676
## 127  0.906462314   -1.1283948
## 128  0.755595077   -1.0072676
## 134  1.057329550   -0.7650132
## 139  0.604727841   -1.1283948
## 142  1.962532967   -0.7650132
## 146  1.660798494   -0.6438861
## 147  1.057329550   -0.8861404

# Tabela de confusão.
table(m1@ymatrix, m1@fitted)

##     
##      -1  1
##   -1 46  4
##   1   3 47

# Classificações incorretas.
i <- which(m1@ymatrix * m1@fitted < 0)

# Fazendo a predição.
grid$y_ksvm <- predict(m1, 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 = c("orange", "cyan")[as.integer(grid$y_svm)],
       pch = 3,
       cex = 0.5)
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.6,
           col = "yellow")
}

3 Funções kernel

3.1 Explorando as funções kernel

Essa seção explora um pouco das funções kernel.

# Dois vetores para empregar nas funções kernel.
x1 <- c(-1, 0, 1)
x2 <- c(0, 1, 2)

# Funções kernel.
grep(ls("package:kernlab"), pattern = "dot$", value = TRUE)

## [1] "anovadot"   "besseldot"  "laplacedot" "polydot"    "rbfdot"    
## [6] "splinedot"  "stringdot"  "tanhdot"    "vanilladot"

# Kernel linear (default)
knl <- vanilladot()
kpar(knl)   # Parâmetros da função kernel.

## list()

knl(x1, x2) # Aplicação.

##      [,1]
## [1,]    2

# Kernel polinomial.
knl <- polydot(degree = 3, scale = 1, offset = 1)
kpar(knl)   # Parâmetros da função kernel.

## $degree
## [1] 3
## 
## $scale
## [1] 1
## 
## $offset
## [1] 1

knl(x1, x2) # Aplicação.

##      [,1]
## [1,]   27

# Kernel de base radial.
knl <- rbfdot(sigma = 1)
kpar(knl)   # Parâmetros da função kernel.

## $sigma
## [1] 1

knl(x1, x2) # Aplicação.

##            [,1]
## [1,] 0.04978707

# Para exibir o código da função kernel.
str(knl)

## Formal class 'rbfkernel' [package "kernlab"] with 2 slots
##   ..@ .Data:function (x, y = NULL)  
##   ..@ kpar :List of 1
##   .. ..$ sigma: num 1

# knl@.Data
body(knl)

## {
##     if (!is(x, "vector")) 
##         stop("x must be a vector")
##     if (!is(y, "vector") && !is.null(y)) 
##         stop("y must a vector")
##     if (is(x, "vector") && is.null(y)) {
##         return(1)
##     }
##     if (is(x, "vector") && is(y, "vector")) {
##         if (!length(x) == length(y)) 
##             stop("number of dimension must be the same on both data points")
##         return(exp(sigma * (2 * crossprod(x, y) - crossprod(x) - 
##             crossprod(y))))
##     }
## }

class(knl)

## [1] "rbfkernel"
## attr(,"package")
## [1] "kernlab"

# Função kernel definida pelo usuário (tirado da documentação).
knl <- function(x, y) {
    (sum(x * y) + 1) * exp(-0.001 * sum((x - y)^2))
}
class(knl) <- "kernel"

# Avaliando a função kernel.
knl(x1, x2)

## [1] 2.991013

#-----------------------------------------------------------------------
# Expansão de características por meio de funções kernel.

# help(kernelMatrix, h = "html")

# Duas características.
x1 <- cbind(head(iris$Sepal.Length)) # Vetor coluna.
x2 <- cbind(head(iris$Petal.Length)) # Vetor coluna.
cbind(x1, x2)

##      [,1] [,2]
## [1,]  5.1  1.4
## [2,]  4.9  1.4
## [3,]  4.7  1.3
## [4,]  4.6  1.5
## [5,]  5.0  1.4
## [6,]  5.4  1.7

# Kernel linear.
# body(vanilladot())
kernelMatrix(kernel = vanilladot(), x1, x2)

## An object of class "kernelMatrix"
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 7.14 7.14 6.63 7.65 7.14 8.67
## [2,] 6.86 6.86 6.37 7.35 6.86 8.33
## [3,] 6.58 6.58 6.11 7.05 6.58 7.99
## [4,] 6.44 6.44 5.98 6.90 6.44 7.82
## [5,] 7.00 7.00 6.50 7.50 7.00 8.50
## [6,] 7.56 7.56 7.02 8.10 7.56 9.18

x1 %*% t(x2)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 7.14 7.14 6.63 7.65 7.14 8.67
## [2,] 6.86 6.86 6.37 7.35 6.86 8.33
## [3,] 6.58 6.58 6.11 7.05 6.58 7.99
## [4,] 6.44 6.44 5.98 6.90 6.44 7.82
## [5,] 7.00 7.00 6.50 7.50 7.00 8.50
## [6,] 7.56 7.56 7.02 8.10 7.56 9.18

tcrossprod(x1, x2)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 7.14 7.14 6.63 7.65 7.14 8.67
## [2,] 6.86 6.86 6.37 7.35 6.86 8.33
## [3,] 6.58 6.58 6.11 7.05 6.58 7.99
## [4,] 6.44 6.44 5.98 6.90 6.44 7.82
## [5,] 7.00 7.00 6.50 7.50 7.00 8.50
## [6,] 7.56 7.56 7.02 8.10 7.56 9.18

# Kernel polinomial.
body(polydot(degree = 1, scale = 1, offset = 1))

## {
##     if (!is(x, "vector")) 
##         stop("x must be a vector")
##     if (!is(y, "vector") && !is.null(y)) 
##         stop("y must be a vector")
##     if (is(x, "vector") && is.null(y)) {
##         (scale * crossprod(x) + offset)^degree
##     }
##     if (is(x, "vector") && is(y, "vector")) {
##         if (!length(x) == length(y)) 
##             stop("number of dimension must be the same on both data points")
##         (scale * crossprod(x, y) + offset)^degree
##     }
## }

# Radial base function kernel.
body(rbfdot(sigma = 2))

## {
##     if (!is(x, "vector")) 
##         stop("x must be a vector")
##     if (!is(y, "vector") && !is.null(y)) 
##         stop("y must a vector")
##     if (is(x, "vector") && is.null(y)) {
##         return(1)
##     }
##     if (is(x, "vector") && is(y, "vector")) {
##         if (!length(x) == length(y)) 
##             stop("number of dimension must be the same on both data points")
##         return(exp(sigma * (2 * crossprod(x, y) - crossprod(x) - 
##             crossprod(y))))
##     }
## }

Dê uma olhada nestes materiais: https://cran.r-project.org/web/packages/kernlab/vignettes/kernlab.pdf, https://quantsignals.wordpress.com/2012/09/25/kernel-learning-svm/, http://www.eric-kim.net/eric-kim-net/posts/1/kernel_trick.html.

# Para abrir o código fonte da função (classe S4).
ksvm
showMethods(ksvm)

getMethod(f = ksvm, signature = "formula")
getMethod(f = ksvm, signature = "matrix")

3.2 Ajustes com diferentes funções kernel

3.2.1 Linear

# k <- grep(ls("package:kernlab"), pattern = "dot$", value = TRUE)
# dput(k)
k <- c("besseldot",
       "laplacedot",
       "polydot",
       "rbfdot",
       "splinedot",
       "tanhdot",
       "vanilladot")

# Todos os ajustes com o valor default dos parâmetros da função kernel.
fits <- sapply(k,
               simplify = FALSE,
               FUN = ksvm,
               x = X,
               y = y,
               scale = FALSE,
               type = "C-svc",
               C = 1)

##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters  
##  Setting default kernel parameters

# Valor da função objetivo (ordenado) e núm de pontos de suporte.
perf <- data.frame(obj = sapply(fits, FUN = getElement, "obj"),
                   nSV = sapply(fits, FUN = getElement, "nSV"))
perf[order(-perf$obj), ]

##                   obj nSV
## laplacedot  -19.64076  55
## rbfdot      -20.64337  39
## polydot     -20.85885  26
## vanilladot  -20.85885  26
## besseldot   -31.58996  42
## tanhdot     -91.17053  50
## splinedot  -949.56121  79

# Acurácia.
perf$acc <- sapply(fits,
                   FUN = function(m0) {
                       # Tabela de confusão.
                       tb <- table(m0@ymatrix, m0@fitted)
                       sum(diag(tb))/sum(tb)
                   })
perf[order(-perf$obj), ]

##                   obj nSV  acc
## laplacedot  -19.64076  55 0.98
## rbfdot      -20.64337  39 0.95
## polydot     -20.85885  26 0.93
## vanilladot  -20.85885  26 0.93
## besseldot   -31.58996  42 0.94
## tanhdot     -91.17053  50 0.68
## splinedot  -949.56121  79 0.25

# Extração dos valores preditos.
fitteds <- lapply(fits, FUN = predict, newdata = grid[, 1:2])
str(fitteds)

## List of 7
##  $ besseldot : num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ laplacedot: num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ polydot   : num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ rbfdot    : num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ splinedot : num [1:3721] 1 1 1 1 1 1 1 1 1 1 ...
##  $ tanhdot   : num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ vanilladot: num [1:3721] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...

# Cria tabela com eles.
fitteds <- data.frame(kernel = rep(k, each = nrow(grid)),
                      pred = unlist(fitteds))
fitteds <- cbind(fitteds, grid[, 1:2])

# Ordena as funções kernel pela acurácia.
fitteds$kernel <- factor(fitteds$kernel,
                         levels = rownames(perf)[order(-perf$acc)])
levels(fitteds$kernel)

## [1] "laplacedot" "rbfdot"     "besseldot"  "polydot"    "vanilladot"
## [6] "tanhdot"    "splinedot"

str(fitteds)

## 'data.frame':    26047 obs. of  4 variables:
##  $ kernel      : Factor w/ 7 levels "laplacedot","rbfdot",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ pred        : num  -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
##  $ Sepal.Length: num  -1.055 -0.979 -0.904 -0.829 -0.753 ...
##  $ Petal.Length: num  -3.31 -3.31 -3.31 -3.31 -3.31 ...

# Gráfico.
xyplot(Petal.Length ~ Sepal.Length | kernel,
       groups = pred,
       data = fitteds,
       as.table = TRUE,
       aspect = 1,
       pch = 1,
       alpha = 0.5,
       cex = 0.5) +
    as.layer({
        xyplot(X[, 2] ~ X[, 1],
               groups = y,
               pch = 19)
    })

Na prática, os hiperparâmetros das funções são determinados por valiação cruzada, bem como o parâmetro de penalidade C.

Neste tutorial foi considerado o SVM apenas para classificação, mas versões para regressão e descoberta de novidade também estão disponíveis. Consulte a documentação e veja o argumento type da ksmv().

Machine Learning

`leg.ufpr.br/~walmes/ensino/ML`

Máquinas de vetores de suporte

Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira

2018-10-17

1 Classificador de vetores de suporte

1.1 Simulando um caso linearmente separável

1.2 Transformando em não linearmente separável

1.3 Otimização quadrática com restrições

2 Implementações de máquinas de vetores de suporte

2.1 Pacote `e1071`

2.2 Pacote `kernlab`

3 Funções kernel

3.1 Explorando as funções kernel

3.2 Ajustes com diferentes funções kernel

3.2.1 Linear

Machine Learning

leg.ufpr.br/~walmes/ensino/ML

Máquinas de vetores de suporte

Prof. Walmes M. Zeviani & Prof. Eduardo V. Ferreira

2018-10-17

1 Classificador de vetores de suporte

1.1 Simulando um caso linearmente separável

1.2 Transformando em não linearmente separável

1.3 Otimização quadrática com restrições

2 Implementações de máquinas de vetores de suporte

2.1 Pacote e1071

2.2 Pacote kernlab

3 Funções kernel

3.1 Explorando as funções kernel

3.2 Ajustes com diferentes funções kernel

3.2.1 Linear

`leg.ufpr.br/~walmes/ensino/ML`

2.1 Pacote `e1071`

2.2 Pacote `kernlab`