Máquinas de vectores de soporte

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs

# 40 puntos sintéticos en dos grupos, 20 de cada clase
X, y = make_blobs(n_samples=40, centers=[[-1.5, -1], [1.5, 1]],
                  cluster_std=0.8, random_state=1)
# Las SVM usan las clases -1 y +1 en lugar de 0 y 1
signo = np.where(y == 1, 1, -1)
print("Clases como 0 y 1:  ", y[:8])
print("Clases como signos: ", signo[:8])

def dibujar_clases(X, y, ax=None, s=30):
    """Clase 0: círculos de color C0. Clase 1: cuadrados de color C1."""
    ax = ax or plt.gca()
    ax.scatter(X[y == 0, 0], X[y == 0, 1], s=s, color="C0",
               marker="o", label="Clase 0")
    ax.scatter(X[y == 1, 0], X[y == 1, 1], s=s, color="C1",
               marker="s", label="Clase 1")

# Tres rectas w1·x1 + w2·x2 + b = 0: nombre, (w1, w2), b
# y dónde escribir el nombre en el dibujo
rectas = [("A", (1, 1), 0, (-2.9, 3)),
          ("B", (3, 1), 0, (-0.8, 3)),
          ("C", (1, 3), 0.5, (3, -0.9))]
x1 = np.array([-4, 4])
dibujar_clases(X, y)
for nombre, (w1, w2), b, (tx, ty) in rectas:
    # Despejamos x2 de la ecuación de la recta
    plt.plot(x1, -(w1 * x1 + b) / w2, color="C2")
    plt.text(tx, ty, nombre, color="C2", fontsize=16)
plt.xlim(-3.5, 3.5)
plt.ylim(-3.5, 3.5)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("Tres rectas que separan las clases")
plt.show()

def distancias(X, w, b):
    """Distancia de cada fila de X a la recta w·x + b = 0."""
    return np.abs(X @ w + b) / np.linalg.norm(w)

for nombre, (w1, w2), b, _ in rectas:
    w = np.array([w1, w2])
    aciertos = np.mean(np.sign(X @ w + b) == signo)
    cerca = distancias(X, w, b).min()
    print(f"Recta {nombre}, w = ({w1}, {w2}), b = {b}:  aciertos"
          f" {aciertos:.0%}   punto más cercano a {cerca:.2f}")

yz = np.linspace(-2, 3, 300)          # la clase por el puntaje
bisagra = np.maximum(0, 1 - yz)
entropia = np.log(1 + np.exp(-yz))

for valor in [-1, 0, 0.5, 1, 2, 3]:
    print(f"y·z = {valor:4}:  bisagra {max(0, 1 - valor):.2f}"
          f"   entropía cruzada {np.log(1 + np.exp(-valor)):.2f}")

plt.plot(yz, bisagra, color="C1", linewidth=2.6, label="Bisagra")
plt.plot(yz, entropia, color="C2", linestyle="--",
         label="Entropía cruzada")
plt.axvline(1, color="C2", linewidth=1)
plt.xlabel("y·z  (positivo: lado correcto)")
plt.ylabel("Pérdida")
plt.title("La bisagra no cobra fuera de la calle")
# Con fondo blanco, la línea vertical no cruza el texto de la leyenda
plt.legend(frameon=True)
plt.show()

def objetivo(X, signo, w, b, C):
    """J = ½‖w‖² + C · (suma de las pérdidas bisagra)."""
    multas = np.maximum(0, 1 - signo * (X @ w + b))
    return 0.5 * w @ w + C * multas.sum()

def svm_subgradiente(X, signo, C, tasa, pasos):
    """SVM lineal con descenso de subgradiente.

    Devuelve los pesos w, el término independiente b y el objetivo J
    antes de cada paso.
    """
    n, d = X.shape
    w = np.zeros(d)                   # un peso por característica
    b = 0.0
    historial = []
    for paso in range(pasos):
        historial.append(objetivo(X, signo, w, b, C))
        yz = signo * (X @ w + b)      # la clase por el puntaje
        multa = yz < 1                # los que se meten en la calle
        # X[multa].T @ signo[multa] es la suma de signo[i] * X[i]
        # sobre los puntos que pagan multa (.T es la transpuesta)
        dw = w - C * X[multa].T @ signo[multa]
        db = -C * signo[multa].sum()
        w = w - tasa * dw             # un paso en contra
        b = b - tasa * db
    return w, b, np.array(historial)

w_svm, b_svm, historial = svm_subgradiente(X, signo, C=1, tasa=0.001,
                                           pasos=20000)
for paso in [0, 10, 100, 1000, 19999]:
    print(f"Paso {paso:5d}:   J = {historial[paso]:.4f}")
print("w =", np.round(w_svm, 3), "  b =", round(b_svm, 3))

from sklearn.svm import SVC

svm = SVC(kernel="linear", C=1).fit(X, y)
print("scikit-learn:   w =", np.round(svm.coef_[0], 3),
      "  b =", np.round(svm.intercept_[0], 3))
print(f"Ancho de la calle: {2 / np.linalg.norm(w_svm):.2f}")

# Vectores de soporte: los que tocan la calle o se meten en ella
yz = signo * (X @ w_svm + b_svm)
nuestros = np.flatnonzero(yz < 1.01)
print("Vectores de soporte (nuestro):     ", nuestros)
print("Vectores de soporte (scikit-learn):", np.sort(svm.support_))
print("Su y·z:", np.round(yz[nuestros], 2))

def dibujar_svm(w, b, ax=None, vectores=None):
    """La frontera (z = 0), los bordes (z = -1 y z = 1) y la calle."""
    ax = ax or plt.gca()
    A, B = np.meshgrid(np.linspace(-4, 4, 300), np.linspace(-4, 4, 300))
    Z = w[0] * A + w[1] * B + b
    ax.contourf(A, B, Z, levels=[-1, 1], colors=["C2"], alpha=0.18)
    ax.contour(A, B, Z, levels=[-1, 0, 1], colors=["C2", "C0", "C2"],
               linestyles=["--", "-", "--"], linewidths=[1.4, 2.6, 1.4])
    if vectores is not None:
        # Un anillo alrededor de cada vector de soporte
        ax.scatter(vectores[:, 0], vectores[:, 1], s=190,
                   facecolors="none", edgecolors="C0", linewidths=1.5)
    ax.set_xlim(-3.5, 3.5)
    ax.set_ylim(-3.5, 3.5)

dibujar_svm(w_svm, b_svm, vectores=X[nuestros])
dibujar_clases(X, y)
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
plt.title("La calle más ancha (C = 1)")
plt.show()

# Un círculo de la clase 0 metido en la zona de la clase 1
X_atip = np.vstack([X, [[0.3, 1.0]]])
y_atip = np.append(y, 0)

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for eje, C in zip(ejes.flat, [0.01, 0.1, 1, 100]):
    modelo = SVC(kernel="linear", C=C).fit(X_atip, y_atip)
    w, b = modelo.coef_[0], modelo.intercept_[0]
    ancho = 2 / np.linalg.norm(w)
    dibujar_svm(w, b, ax=eje, vectores=modelo.support_vectors_)
    dibujar_clases(X_atip, y_atip, ax=eje, s=10)
    eje.set_title(f"C = {C}\nancho {ancho:.2f}")
    print(f"C = {C:5}:  ancho {ancho:.2f}"
          f"   vectores de soporte {len(modelo.support_):2d}"
          f"   el atípico va a la clase {modelo.predict([[0.3, 1.0]])[0]}")
plt.show()

from sklearn.datasets import load_wine
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

X_vino, y_vino = load_wine(return_X_y=True)
# Columna 12: prolina. Columna 10: tono del color
print(f"La prolina varía en {np.ptp(X_vino[:, 12]):.0f};"
      f" el tono, en {np.ptp(X_vino[:, 10]):.2f}")

sin_escalar = cross_val_score(SVC(), X_vino, y_vino, cv=5).mean()
escalado = cross_val_score(make_pipeline(StandardScaler(), SVC()),
                           X_vino, y_vino, cv=5).mean()
print(f"Exactitud sin escalar: {sin_escalar:.3f}")
print(f"Exactitud estandarizando: {escalado:.3f}")

from sklearn.datasets import make_circles
from sklearn.svm import LinearSVC

# 200 puntos sintéticos: la clase 1 es el anillo de adentro
X_circ, y_circ = make_circles(n_samples=200, factor=0.4, noise=0.08,
                              random_state=0)
# Levantamos cada punto: la altura es su distancia al cuadrado al centro
x3 = X_circ[:, 0] ** 2 + X_circ[:, 1] ** 2
X_3d = np.column_stack([X_circ, x3])

plano = LinearSVC().fit(X_circ, y_circ)
espacio = LinearSVC().fit(X_3d, y_circ)
print(f"Recta en el plano:     aciertos {plano.score(X_circ, y_circ):.2f}")
print(f"Plano en 3D:           aciertos {espacio.score(X_3d, y_circ):.2f}")
print(f"Altura de la clase 1:  hasta {x3[y_circ == 1].max():.2f}")
print(f"Altura de la clase 0:  desde {x3[y_circ == 0].min():.2f}")

fig = plt.figure()
eje = fig.add_subplot(projection="3d")
for clase, color, marca in [(0, "C0", "o"), (1, "C1", "s")]:
    m = y_circ == clase
    eje.scatter(X_circ[m, 0], X_circ[m, 1], x3[m], color=color,
                marker=marca, s=12, depthshade=False)
# Un plano horizontal a media altura separa las dos clases
A, B = np.meshgrid(np.linspace(-1.2, 1.2, 2), np.linspace(-1.2, 1.2, 2))
eje.plot_surface(A, B, np.full_like(A, 0.48), color="C2", alpha=0.3)
eje.set_xlabel("$x_1$")
eje.set_ylabel("$x_2$")
eje.set_zlabel("$x_3$")
eje.set_title("Levantar los anillos a 3D")
plt.show()

def levantar(p):
    """Tres características nuevas a partir de (x1, x2)."""
    return np.array([p[0] ** 2, np.sqrt(2) * p[0] * p[1], p[1] ** 2])

a = np.array([1.0, 2.0])
c = np.array([3.0, -1.0])
print("Levantando y multiplicando:", round(levantar(a) @ levantar(c), 6))
print("Con el kernel (a·c)²:      ", (a @ c) ** 2)

for nombre, modelo in [("lineal", SVC(kernel="linear")),
                       ("polinomial, grado 2", SVC(kernel="poly",
                                                   degree=2)),
                       ("RBF", SVC(kernel="rbf"))]:
    modelo.fit(X_circ, y_circ)
    print(f"Kernel {nombre:20s} aciertos {modelo.score(X_circ, y_circ):.2f}")

from sklearn.datasets import make_moons
from sklearn.model_selection import train_test_split

# Las dos lunas del tema 8, con la misma división
X_lunas, y_lunas = make_moons(n_samples=200, noise=0.3, random_state=0)
X_entr, X_val, y_entr, y_val = train_test_split(
    X_lunas, y_lunas, test_size=0.3, random_state=0)

A, B = np.meshgrid(np.linspace(-2, 3, 200), np.linspace(-1.5, 2, 200))
rejilla = np.column_stack([A.ravel(), B.ravel()])

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for eje, (gamma, C) in zip(ejes.flat,
                           [(0.1, 1), (0.1, 100), (10, 1), (10, 100)]):
    modelo = SVC(kernel="rbf", C=C, gamma=gamma).fit(X_entr, y_entr)
    zona = modelo.predict(rejilla).reshape(A.shape)
    # Cada región, del color de su clase y muy clara
    eje.contourf(A, B, zona, levels=[-0.5, 0.5, 1.5],
                 colors=["C0", "C1"], alpha=0.2)
    dibujar_clases(X_entr, y_entr, ax=eje, s=8)
    eje.set_title(f"gamma = {gamma}\nC = {C}")
    print(f"γ = {gamma:4}, C = {C:3}:  aciertos al entrenar "
          f"{modelo.score(X_entr, y_entr):.2f}"
          f"   en validación {modelo.score(X_val, y_val):.2f}")
plt.show()

from sklearn.model_selection import GridSearchCV
import pandas as pd

pipe = make_pipeline(StandardScaler(), SVC(kernel="rbf"))
rejilla_hp = {"svc__C": [0.1, 1, 10, 100],
              "svc__gamma": [0.1, 1, 10, 100]}
busqueda = GridSearchCV(pipe, rejilla_hp, cv=5)
busqueda.fit(X_entr, y_entr)

# La exactitud media de validación cruzada de las 16 combinaciones
tabla = pd.DataFrame(busqueda.cv_results_["params"])
tabla["exactitud"] = busqueda.cv_results_["mean_test_score"]
print(tabla.pivot(index="svc__C", columns="svc__gamma",
                  values="exactitud").round(2))
print("Mejor combinación:", busqueda.best_params_)
print(f"Exactitud en validación cruzada: {busqueda.best_score_:.3f}")
print(f"Exactitud en el conjunto de validación: "
      f"{busqueda.score(X_val, y_val):.3f}")

import time
from sklearn.datasets import make_classification

X_grande, y_grande = make_classification(n_samples=5000, n_features=20,
                                         random_state=0)
inicio = time.perf_counter()
lento = SVC(kernel="linear").fit(X_grande, y_grande)
t_svc = time.perf_counter() - inicio

inicio = time.perf_counter()
rapido = LinearSVC().fit(X_grande, y_grande)
t_lineal = time.perf_counter() - inicio

print(f"SVC:       aciertos {lento.score(X_grande, y_grande):.3f}"
      f"   vectores de soporte {len(lento.support_)}")
print(f"LinearSVC: aciertos {rapido.score(X_grande, y_grande):.3f}")
print("¿LinearSVC fue más de 10 veces más rápido?", t_svc > 10 * t_lineal)