Support vector machines

Laboratorios

En R, podemos aplicar los conceptos aprendidos mediante la función svm de la librería e1071. Su uso básico es de la forma:

svm(formula, data, kernel, cost)

donde:

• kernel: tipo de núcleo a utilizar.

  • 'linear': \(K(x,x')=\langle x,x'\rangle\)

  • 'polynomial': \(K(x,x') = \left(\gamma\langle x,x'\rangle+b\right)^d\)

  • 'radial': \(K(x,x')=\exp\{-\gamma\|x-x'\|^2\}\)

  • 'sigmoid': \(K(x,x')=\tanh\left(\gamma\langle x,x'\rangle+b\right)\)

• cost: constante \(C\) de regularización.

Por defecto, kernel = 'radial', gamma = 1/(data dimension), cost = 1.

Ejemplos

Hiperplano separador óptimo

Generamos datos linealmente separables:

set.seed(1)

x = matrix(rnorm(40), ncol = 2)
y = c(rep(-1,10), rep(1,10))
x[y==1,] = x[y==1,] + 2

datosA = data.frame(x = x, y = as.factor(y))


# Plot:
figureA1 = ggplot(datosA) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$'), fill = NULL) +
  scale_color_brewer(palette = 'Set1')

Obtenemos el hiperplano separador óptimo aplicando svm con núcleo lineal y \(C\) lo suficientemente grande para no permitir puntos del lado incorrecto del margen.

modeloA = svm(y~., datosA, kernel = 'linear', cost = 1e5, scale = FALSE)

• Resumen del modelo:

summary(modeloA)

## 
## Call:
## svm(formula = y ~ ., data = datosA, kernel = "linear", cost = 1e+05, 
##     scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1e+05 
## 
## Number of Support Vectors:  3
## 
##  ( 1 2 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1

• Identificación de vectores soporte:

modeloA$index

## [1]  2 14 15

• Gráfica:

# Vectores soporte:
supportA = datosA[modeloA$index,]


# Grilla (predicción):
x1grid = seq(min(datosA$x.1), max(datosA$x.1), by = 0.01)
x2grid = seq(min(datosA$x.2), max(datosA$x.2), by = 0.01)
aux = as.matrix(expand.grid(x1grid, x2grid))
gridA = data.frame(x.1 = aux[,1], x.2 = aux[,2])

values = predict(modeloA, gridA, decision.values = TRUE)
values = attr(values, 'decision.values')

gridA$z = values


# Plot:
figureA2 = ggplot(datosA) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, alpha = 0.4, show.legend = F) + 
  geom_contour(data = gridA, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  geom_contour(data = gridA, aes(x = x.1, y = x.2, z = z), breaks = c(-1,1), col = 'gray25', size = 0.5, linetype = 'dashed') +
  geom_point(data = supportA, aes(x = x.1, y = x.2, fill = y), shape = 21, stroke = 1.2, size = 3, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

Clasificador de vectores soporte

Generamos datos linealmente no separables:

set.seed(1)

x = matrix(rnorm(80), ncol = 2)
y = c(rep(-1,20), rep(1,20))
x[y==1,] = x[y==1,] + 2

datosB = data.frame(x = x, y = as.factor(y))


# Plot:
figureB = ggplot(datosB) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1')

Obtenemos el clasificador de vectores soporte aplicando svm con núcleo lineal e incluyendo un costo \(C\): en primer lugar utilizamos el valor por defecto \(C=1\).

modeloB = svm(y~., datosB, kernel = 'linear', scale = FALSE)

• Resumen del modelo:

summary(modeloB)

## 
## Call:
## svm(formula = y ~ ., data = datosB, kernel = "linear", scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
## 
## Number of Support Vectors:  13
## 
##  ( 6 7 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  -1 1

• Identificación de vectores soporte:

modeloB$index

##  [1]  4  8 11 15 16 19 24 27 28 32 34 35 37

• Gráfica:

# Vectores soporte:
supportB = datosB[modeloB$index,]


# Grilla (Predicción):
x1grid = seq(min(datosB$x.1), max(datosB$x.1), by = 0.01)
x2grid = seq(min(datosB$x.2), max(datosB$x.2), by = 0.01)
aux = as.matrix(expand.grid(x1grid, x2grid))
gridB = data.frame(x.1 = aux[,1], x.2 = aux[,2])

values = predict(modeloB, gridB, decision.values = TRUE)
values = attr(values, 'decision.values')

gridB$z = values


# Plot:
figureB2 = ggplot(datosB) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, alpha = 0.4, show.legend = F) + 
  geom_contour(data = gridB, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  geom_contour(data = gridB, aes(x = x.1, y = x.2, z = z), breaks = c(-1,1), col = 'gray25', size = 0.5, linetype = 'dashed') +
  geom_point(data = supportB, aes(x = x.1, y = x.2, fill = y), shape = 21, stroke = 1.2, size = 3, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

---

Efectos del valor de \(C\)

---

Veremos a continuación como afecta a la cantidad de vectores soporte el valor de \(C\).

C = c(0.01, 0.1, 1, 10)
plots = list()

for (i in 1:length(C)){
  modelo = svm(y~., datosB, kernel ='linear', cost = C[i], scale = FALSE)
  
  # Vectores soporte:
  support = datosB[modelo$index,]
  
  # Grilla (predicción):
  grid = gridB
  values = predict(modelo, gridB, decision.values = TRUE)
  values = attr(values,'decision.values')
  grid$z = values
  
  # Plot:
  figure = ggplot(datosB) + 
    geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, alpha = 0.4, show.legend = F) + 
  geom_contour(data = grid, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  geom_contour(data = grid, aes(x = x.1, y = x.2, z = z), breaks = c(-1,1), col = 'gray25', size = 0.5, linetype = 'dashed') +
  geom_point(data = support, aes(x = x.1, y = x.2, fill = y), shape = 21, stroke = 1.2, size = 3, show.legend = F) +
  theme_classic() +
  labs(title = paste('C =',C[i]), subtitle = paste('Vectores soporte:',length(modelo$index)), x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')
  
  plots[[i]] = figure
}

Máquina de vectores soporte

En los dos casos anteriores, hemos visto como utilizar la función svm para hallar fronteras de clasificación lineales. También hemos analizado los efectos del valor de \(C\) respecto a los vectores soporte. Pero, ¿cómo hacemos para elegir un valor de \(C\) adecuado?

Los parámetros en juego necesitan ser ajustados mediante algun proceso de validación. Además, es importante analizar la performance de nuestro clasificador en un nuevo conjunto de datos.

Por lo tanto, ahora que analizaremos el caso general para hallar fronteras de clasificación no lineales, utilizaremos los siguientes recursos:

• Un conjunto de datos de entrenamiento.
• Un proceso de validación para elegir los parámetros con los cuales entrenaremos nuestro clasificador.
• Un conjunto de datos de prueba para evaluar dicho clasificador.

Generamos los datos y los separamos en entrenamiento-prueba:

set.seed(1)

x = matrix(rnorm (400), ncol=2)
x[1:100,] = x[1:100,]+2
x[101:150,] = x[101:150,]-2
y = c(rep(1,150) ,rep(2,50))

datosC = data.frame(x = x, y = as.factor(y))

train = sample(200,120)
datosCtrain = datosC[train,]
datosCtest = datosC[-train,]


# Plot:
figureC1 = ggplot(datosCtrain) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) +
  labs(title = 'Datos de entrenamiento', x = TeX('$x_1$'), y = TeX('$x_2$'), fill = NULL) +
  coord_cartesian(xlim = c(min(datosC$x.1), max(datosC$x.1)), ylim = c(min(datosC$x.2), max(datosC$x.2))) +
  theme_classic() +
  scale_color_brewer(palette = 'Set1')

figureC2 = ggplot(datosCtest) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) +
  labs(title = 'Datos de prueba', x = TeX('$x_1$'), y = TeX('$x_2$'), fill = NULL) +
  coord_cartesian(xlim = c(min(datosC$x.1), max(datosC$x.1)), ylim = c(min(datosC$x.2), max(datosC$x.2))) +
  theme_classic() +
  scale_color_brewer(palette = 'Set1')

En primer lugar, entrenamos una máquina de vectores soporte con núcleo radial (\(\gamma=1/2\)) y costo \(C=1\), valores por defecto de la función svm.

modeloC1 = svm(y~., datosCtrain, scale = FALSE)

• Resumen del modelo:

summary(modeloC1)

## 
## Call:
## svm(formula = y ~ ., data = datosCtrain, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  43
## 
##  ( 24 19 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  1 2

• Identificación de vectores soporte:

modeloC1$index

##  [1]  11  19  20  23  24  27  28  29  43  47  53  61  67  71  73  74  77  80  85
## [20]  97 111 114 116 120   2  12  17  21  26  31  40  56  57  59  65  70  79  81
## [39]  94 101 104 113 117

• Errores:

# Error en datos de entrenamiento:
Ytrainfit = predict(modeloC1, datosCtrain)
errortrain = mean(datosCtrain$y!=Ytrainfit)


# Error en datos de prueba:
Ytestfit = predict(modeloC1, datosCtest)
errortest = mean(datosCtest$y!=Ytestfit)

## Error de entrenamiento        Error de prueba 
##             0.05833333             0.15000000

• Gráfica:

# Vectores soporte:
supportC = datosCtrain[modeloC1$index,]


# Grilla (predicción):
x1.grid = seq(min(datosC$x.1), max(datosC$x.1), by = 0.01)
x2.grid = seq(min(datosC$x.2), max(datosC$x.2), by = 0.01)
aux = as.matrix(expand.grid(x1.grid, x2.grid))
gridC = data.frame(x.1 = aux[,1], x.2 = aux[,2])

values = predict(modeloC1, gridC, decision.values = TRUE)
values = attr(values,'decision.values')

gridC$z = values


# Plot:
figureC3 = ggplot(datosCtrain) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, alpha = 0.4, show.legend = F) + 
  geom_contour(data = gridC, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  geom_contour(data = gridC, aes(x = x.1, y = x.2, z = z), breaks = c(-1,1), col = 'gray25', size = 0.5, linetype = 'dashed') +
  geom_point(data = supportC, aes(x = x.1, y = x.2, fill = y), shape = 21, stroke = 1.2, size = 3, show.legend = F) +
  theme_classic() +
  labs(title = 'Datos de entrenamiento', subtitle = paste('Error :', round(errores[1,1],4)), x = TeX('$x_1$'), y = TeX('$x_2$')) +
  coord_cartesian(xlim = c(min(datosC$x.1), max(datosC$x.1)), ylim = c(min(datosC$x.2), max(datosC$x.2))) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

figureC4 = ggplot(datosCtest) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) + 
  geom_contour(data = gridC, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  theme_classic() +
  labs(title = 'Datos de prueba', subtitle = paste('Error:', round(errores[1,2],4)), x = TeX('$x_1$'), y = TeX('$x_2$')) +
  coord_cartesian(xlim = c(min(datosC$x.1), max(datosC$x.1)), ylim = c(min(datosC$x.2), max(datosC$x.2))) +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

---

Efectos del valor de \(C\)

---

Veremos a continuación como afecta el valor de \(C\).

---

Efectos del valor de \(\gamma\)

---

Veremos a continuación como afecta el valor de \(\gamma\).

---

Elección de valores óptimos para los parámetros

---

A partir de una lista de valores posibles para los parámetros, podemos hacer nuestra elección a partir de un proceso de validación cruzada mediante la función tune del paquete e1071.

set.seed(1)

tuneCV = tune(svm , y∼., data = datosCtrain, ranges = list(cost=c(0.1,1,10,100,1000), gamma=c(0.5,1,2,3,4)))

• Mejores parámetros:

tuneCV$best.parameters

	cost	gamma
2	1	0.5

Finalizamos entonces entrenando una máquina de vectores soporte con núcleo radial (\(\gamma=0.5\)) y \(C=1\), y reportamos el resultado para los datos de prueba.

modeloC2 = svm(y~., datosCtrain, gamma = 0.5, scale = FALSE)


#Error en datos de prueba:
Ytestfit = predict(modeloC2, datosCtest)
errortest = mean(datosCtest$y!=Ytestfit)


# Grilla (predicción):
values = predict(modeloC2, gridC, decision.values = TRUE)
values = attr(values,'decision.values')
gridC$z = values


# Plot:
figureC5 = ggplot(datosCtest) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) + 
  geom_contour(data = gridC, aes(x = x.1, y = x.2, z = z), breaks = 0, col = 'black', size = 1) +
  theme_classic() +
  labs(title = 'Datos de prueba', subtitle = paste('Error:', round(errortest,4)), x = TeX('$x_1$'), y = TeX('$x_2$')) +
  coord_cartesian(xlim = c(min(datosC$x.1), max(datosC$x.1)), ylim = c(min(datosC$x.2), max(datosC$x.2))) +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

---

Curva ROC

---

Para dibujar curvas ROC se puede utilizar la librería ROCR. Trazaremos dichas curvas para el modelo óptimo, tanto para los datos de entrenamiento como para los datos de prueba.

# ROC train:
Ytrainvalues = predict(modeloC2, datosCtrain, decision.values = TRUE)
Ytrainvalues = attr(Ytrainvalues, 'decision.values')

predtrain = prediction(Ytrainvalues, datosCtrain$y, label.ordering = c('2', '1'))
perftrain = performance(predtrain, 'tpr', 'fpr')
datosROCtrain = data.frame(FPR = perftrain@x.values[[1]], TPR = perftrain@y.values[[1]])

AUCtrain = performance(predtrain, 'auc')
AUCtrain = as.numeric(AUCtrain@y.values)


# ROC test:
Ytestvalues = predict(modeloC2, datosCtest, decision.values = TRUE)
Ytestvalues = attr(Ytestvalues, 'decision.values')

predtest = prediction(Ytestvalues, datosCtest$y, label.ordering = c('2', '1'))
perftest = performance(predtest, 'tpr', 'fpr')
datosROCtest = data.frame(FPR = perftest@x.values[[1]], TPR = perftest@y.values[[1]])

AUCtest = performance(predtest, 'auc')
AUCtest = as.numeric(AUCtest@y.values)


# Plots:
figureC6 = ggplot() + 
  geom_line(data = datosROCtrain, aes(x = FPR, y = TPR), color = 'darkblue', size = 1.5) +
  geom_line(aes(x = c(0,1), y = c(0,1)), color = 'gray', linetype = 'dashed') +
  theme_classic() +
  labs(title = 'Curva ROC - Datos de entrenamiento', subtitle = paste('AUC:', round(AUCtrain,2)))

figureC7 = ggplot() + 
  geom_line(data = datosROCtest, aes(x = FPR, y = TPR), color = 'darkred', size = 1.5) +
  geom_line(aes(x = c(0,1), y = c(0,1)), color = 'gray', linetype = 'dashed') +
  theme_classic() +
  labs(title = 'Curva ROC - Datos de prueba', subtitle = paste('AUC:', round(AUCtest,2)))

Finalizamos trazando curvas ROC para diferentes valores del parámetro \(\gamma\) sobre los datos de prueba.

Máquina de vectores soporte multiclase

Para el caso multiclase, la función svm aborda el problema realizando los clasificadores 1-vs-1.

set.seed(1)

x = matrix(rnorm (500), ncol=2)
x[1:100,] = x[1:100,]+2
x[101:150,] = x[101:150,]-2
x[201:250,2] = x[201:250,2]+2
y = c(rep(1,150),rep(2,50), rep(3,50))

datosD = data.frame(x = x, y = as.factor(y))


# Plot:
figureD = ggplot(datosD) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$'), fill = NULL) +
  scale_color_brewer(palette = 'Set1')

Entrenamos una máquina de vectores soporte con núcleo radial (\(\gamma=1/2\) por defecto) y \(C=10\).

modeloD = svm(y~., datosD, cost = 10, scale = FALSE)

• Resumen del modelo:

summary(modeloD)

## 
## Call:
## svm(formula = y ~ ., data = datosD, cost = 10, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  10 
## 
## Number of Support Vectors:  112
## 
##  ( 44 31 37 )
## 
## 
## Number of Classes:  3 
## 
## Levels: 
##  1 2 3

• Identificación de vectores soporte:

modeloD$index

##   [1]   3   7  13  14  20  24  25  27  28  29  35  36  45  52  53  54  56  61
##  [19]  63  65  67  74  75  81  84  89  90  95  97  99 106 107 109 110 111 113
##  [37] 122 130 133 135 141 142 147 149 151 152 155 157 159 160 161 162 163 166
##  [55] 167 169 170 171 174 175 176 178 180 181 182 185 187 189 190 192 194 195
##  [73] 196 199 200 201 202 203 204 205 206 207 208 209 210 211 214 215 216 217
##  [91] 218 219 221 222 225 226 228 229 230 234 236 237 238 239 240 241 242 243
## [109] 245 248 249 250

• Gráfica:

# Vectores soporte:
supportD = datosD[modeloD$index,]


# Grilla:
x1.grid = seq(min(datosD$x.1), max(datosD$x.1), by = 0.1)
x2.grid = seq(min(datosD$x.2), max(datosD$x.2), by = 0.1)
aux = as.matrix(expand.grid(x1.grid, x2.grid))
gridD = data.frame(x.1 = aux[,1], x.2 = aux[,2])

clasif = predict(modeloD, gridD)
gridD$clasif = clasif


# Plot:
figureD = ggplot(datosD) + 
  geom_point(mapping = aes(x = x.1, y = x.2, color = y), shape = 19, size = 3, show.legend = F) + 
  geom_point(data = gridD, aes(x = x.1, y = x.2, color = clasif), shape = 15, size = 2.5, alpha = 0.1, stroke = 0, show.legend = F) +
  theme_classic() +
  labs(x = TeX('$x_1$'), y = TeX('$x_2$')) +
  scale_color_brewer(palette = 'Set1') +
  scale_fill_brewer(palette = 'Set1')

Referencias

James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2013. An Introduction to Statistical Learning. Vol. 112. Springer.