Agrupar sin etiquetas: clustering

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

# 300 puntos sintéticos alrededor de 4 centros
X, y_verdad = make_blobs(n_samples=300, centers=4, cluster_std=0.9,
                         random_state=7)
print("Forma de X:", X.shape)
print("Primeros puntos:", X[:3].round(1).tolist())

plt.scatter(X[:, 0], X[:, 1], s=16)
plt.xlabel("Característica 1")
plt.ylabel("Característica 2")
plt.title("300 puntos sin etiqueta (sintéticos)")
plt.show()

def kmeans(X, centroides, pasos=100):
    """K-means desde unos centroides iniciales.

    Devuelve el grupo de cada punto, los centroides finales
    y la lista de centroides de cada vuelta.
    """
    k = len(centroides)
    historial = [centroides]
    for _ in range(pasos):
        # 1) Asignar: distancia de cada punto a cada centroide, (n, k)
        d = np.linalg.norm(X[:, None, :] - centroides[None, :, :], axis=2)
        grupo = d.argmin(axis=1)
        # 2) Mover: cada centroide va al promedio de su grupo
        #    (si un grupo quedó vacío, su centroide no se mueve)
        nuevos = np.array([X[grupo == j].mean(axis=0) if np.any(grupo == j)
                           else centroides[j] for j in range(k)])
        if np.allclose(nuevos, centroides):   # ya no se mueven: listo
            break
        centroides = nuevos
        historial.append(centroides)
    return grupo, centroides, historial

def inercia(X, grupo, centroides):
    """Suma de las distancias al cuadrado de cada punto a su centroide."""
    return np.sum((X - centroides[grupo]) ** 2)

rng = np.random.default_rng(6)
# Cuatro puntos de los datos, al azar, como centroides iniciales
iniciales = X[rng.choice(len(X), size=4, replace=False)]
grupo, centroides, historial = kmeans(X, iniciales)

print("Vueltas hasta que los centroides se quedaron quietos:",
      len(historial) - 1)
print("Centroides finales:")
print(centroides.round(2))
print(f"Inercia final: {inercia(X, grupo, centroides):.1f}")
print("Puntos por grupo:", np.bincount(grupo))

FORMAS = ["o", "s", "^", "D"]          # círculo, cuadrado, triángulo, rombo
COLORES = ["C0", "C2", "C3", "C0"]

def dibujar_grupos(eje, X, grupo, s=14):
    """Un color y una forma por grupo; el ruido (grupo -1) va como x."""
    for g in np.unique(grupo):
        parte = X[grupo == g]
        if g == -1:
            eje.scatter(parte[:, 0], parte[:, 1], s=s, marker="x",
                        color="C2", linewidths=1)
        else:
            eje.scatter(parte[:, 0], parte[:, 1], s=s, marker=FORMAS[g % 4],
                        color=COLORES[g % 4], edgecolors="C0",
                        linewidths=0.6)

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for vuelta, eje in enumerate(ejes.flat):
    c = historial[vuelta]
    # Asignar cada punto al centroide de esta vuelta
    g = np.linalg.norm(X[:, None, :] - c[None, :, :], axis=2).argmin(axis=1)
    dibujar_grupos(eje, X, g, s=12)
    final = vuelta == len(historial) - 1
    eje.scatter(c[:, 0], c[:, 1], marker="X", s=150, linewidths=1.5,
                color="C1" if final else "C0", edgecolors="white")
    nombre = "Inicio" if vuelta == 0 else f"Vuelta {vuelta}"
    eje.set_title(f"{nombre}\ninercia {inercia(X, g, c):.0f}")
plt.show()

for semilla in range(6):
    rng = np.random.default_rng(semilla)
    inicio = X[rng.choice(len(X), size=4, replace=False)]
    g, c, h = kmeans(X, inicio)
    print(f"Semilla {semilla}:  {len(h) - 1:2d} vueltas   "
          f"inercia {inercia(X, g, c):7.1f}   "
          f"tamaños {np.bincount(g, minlength=4)}")

from sklearn.cluster import KMeans

modelo = KMeans(n_clusters=4, n_init=10, random_state=0)
grupos_km = modelo.fit_predict(X)

print(f"Inercia: {modelo.inertia_:.1f}")
print("Centroides:")
print(modelo.cluster_centers_.round(2))
print("Grupo de los 10 primeros puntos:", grupos_km[:10])

from sklearn.metrics import silhouette_score

ks = range(2, 9)
inercias, siluetas = [], []
for k in ks:
    m = KMeans(n_clusters=k, n_init=10, random_state=0).fit(X)
    inercias.append(m.inertia_)
    siluetas.append(silhouette_score(X, m.labels_))
    print(f"k = {k}:   inercia {m.inertia_:7.1f}   "
          f"silueta {siluetas[-1]:.3f}")

mejor_k = ks[np.argmax(siluetas)]
fig, (arriba, abajo) = plt.subplots(2, 1, figsize=(4.8, 4.8), sharex=True)
arriba.plot(ks, inercias, marker="o", color="C0")
arriba.set_ylabel("Inercia")
arriba.set_title("El codo y la silueta según k")
abajo.plot(ks, siluetas, marker="o", color="C2")
abajo.scatter([mejor_k], [max(siluetas)], s=160, color="C1", zorder=3,
              label=f"La más alta: k = {mejor_k}")
abajo.set_ylabel("Silueta")
abajo.set_xlabel("Cantidad de grupos k")
# Solo para dibujar: margen para que el punto no toque el borde de arriba
abajo.margins(y=0.15)
abajo.legend()
plt.show()

from sklearn.datasets import make_moons
from sklearn.cluster import DBSCAN

X_lunas, y_lunas = make_moons(n_samples=300, noise=0.08, random_state=0)

db = DBSCAN(eps=0.2, min_samples=5).fit(X_lunas)
grupos_db = db.labels_
nucleo = np.zeros(len(X_lunas), dtype=bool)
nucleo[db.core_sample_indices_] = True
ruido = grupos_db == -1
borde = ~nucleo & ~ruido

print("Grupos encontrados:", len(set(grupos_db) - {-1}))
print(f"Núcleo: {nucleo.sum()}   borde: {borde.sum()}   "
      f"ruido: {ruido.sum()}")

from scipy.cluster.hierarchy import linkage, dendrogram

# Cada fila de Z es una unión: grupo, grupo, distancia, tamaño
Z = linkage(X, method="ward")
print("Forma de Z:", Z.shape)
print("Altura de las 4 últimas uniones:", Z[-4:, 2].round(1))

# Un corte entre la 4.ª y la 3.ª unión desde arriba deja 4 grupos
corte = (Z[-4, 2] + Z[-3, 2]) / 2
dendrogram(Z, truncate_mode="lastp", p=12, leaf_font_size=11,
           link_color_func=lambda k: "C0")
plt.axhline(corte, color="C1", linestyle="--", label="Corte: 4 grupos")
plt.ylim(0, 210)                      # espacio arriba para la leyenda
plt.ylabel("Distancia de unión")
plt.title("Dendrograma (enlace ward)")
plt.legend(loc="upper left")
plt.show()

from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics import adjusted_rand_score

grupos_ag = AgglomerativeClustering(n_clusters=4,
                                    linkage="ward").fit_predict(X)
print("Puntos por grupo:", np.bincount(grupos_ag))
print("Acuerdo con K-means:", adjusted_rand_score(grupos_km, grupos_ag))
print("Acuerdo con la verdad:", adjusted_rand_score(y_verdad, grupos_ag))

from sklearn.mixture import GaussianMixture

# return_centers=True devuelve también los centros de las manchas
X_m, y_est, centros = make_blobs(n_samples=300, centers=3,
                                 cluster_std=1.15, random_state=170,
                                 return_centers=True)
estirar = np.array([[0.6, -0.6], [-0.4, 0.8]])
# Cada punto (fila) por la matriz: la nube se estira y se gira
X_est = X_m @ estirar

grupos_km3 = KMeans(n_clusters=3, n_init=10,
                    random_state=0).fit_predict(X_est)
gmm = GaussianMixture(n_components=3, n_init=5, random_state=0)
grupos_gmm = gmm.fit_predict(X_est)

print("Acuerdo con la verdad:")
print(f"  K-means: {adjusted_rand_score(y_est, grupos_km3):.2f}")
print(f"  Mezcla:  {adjusted_rand_score(y_est, grupos_gmm):.2f}")
print("Pesos de las campanas:", gmm.weights_.round(2))

prob = gmm.predict_proba(X_est)
print("Probabilidades de los 3 primeros puntos:")
print(prob[:3].round(3))

seguridad = prob.max(axis=1)        # la probabilidad del grupo elegido
dudoso = seguridad.argmin()
print("El punto más dudoso:", X_est[dudoso].round(2),
      "probabilidades", prob[dudoso].round(2))
print("Puntos con menos de 95 % de seguridad:", (seguridad < 0.95).sum())

# Una rejilla de puntos que cubre los datos, para calcular la densidad
A, B = np.meshgrid(np.linspace(-7, 4, 200), np.linspace(-4, 6, 200))
rejilla = np.column_stack([A.ravel(), B.ravel()])
# score_samples da el logaritmo de la densidad; exp lo deshace
densidad = np.exp(gmm.score_samples(rejilla)).reshape(A.shape)

fig, eje = plt.subplots()
eje.contour(A, B, densidad, levels=6, colors="C2", linewidths=0.8)
dibujar_grupos(eje, X_est, grupos_gmm, s=12)
dudosos = X_est[seguridad < 0.95]
eje.scatter(dudosos[:, 0], dudosos[:, 1], s=140, facecolors="none",
            edgecolors="C1", linewidths=1.8, label="Menos de 95 % seguro")
eje.set_title("Mezcla de 3 gaussianas: elipses")
eje.legend(loc="lower left")
plt.show()

for k in range(1, 7):
    m = GaussianMixture(n_components=k, n_init=5,
                        random_state=0).fit(X_est)
    print(f"k = {k}:   log-verosimilitud media {m.score(X_est):6.3f}"
          f"   BIC {m.bic(X_est):6.1f}")

from sklearn.neighbors import KernelDensity

rng = np.random.default_rng(1)
tiempos = np.concatenate([rng.normal(12, 3, 70), rng.normal(32, 6, 50)])

def densidad_normal(x, mu, sigma):
    """La fórmula de la campana, del tema 5."""
    return (np.exp(-(x - mu) ** 2 / (2 * sigma ** 2))
            / (sigma * np.sqrt(2 * np.pi)))

def kde(xs, datos, h):
    """Una campanita de ancho h sobre cada dato, y el promedio."""
    # Tabla (len(xs), len(datos)); el promedio va por filas
    return densidad_normal(xs[:, None], datos[None, :], h).mean(axis=1)

xs = np.array([5.0, 12.0, 22.0, 32.0])
print("Desde cero:  ", kde(xs, tiempos, h=2).round(4))
# scikit-learn quiere una columna, y da el logaritmo de la densidad
modelo_kde = KernelDensity(bandwidth=2).fit(tiempos.reshape(-1, 1))
log_d = modelo_kde.score_samples(xs.reshape(-1, 1))
print("scikit-learn:", np.exp(log_d).round(4))

from sklearn.model_selection import GridSearchCV, KFold

anchos = [0.3, 1, 2, 3, 4, 6, 10]
busqueda = GridSearchCV(KernelDensity(), {"bandwidth": anchos},
                        cv=KFold(5, shuffle=True, random_state=0))
busqueda.fit(tiempos.reshape(-1, 1))
for h, s in zip(anchos, busqueda.cv_results_["mean_test_score"]):
    print(f"h = {h:4}:  log-verosimilitud en validación {s:7.1f}")
h_mejor = busqueda.best_params_["bandwidth"]
print("Mejor ancho de banda:", h_mejor)

xs = np.linspace(0, 50, 400)
plt.hist(tiempos, bins=25, density=True, color="C3", label="Histograma")
plt.plot(xs, kde(xs, tiempos, 0.3), color="C2", linewidth=1,
         label="h = 0,3")
plt.plot(xs, kde(xs, tiempos, 10), color="C0", linestyle="--",
         label="h = 10")
plt.plot(xs, kde(xs, tiempos, h_mejor), color="C1", linewidth=2.5,
         label=f"h = {h_mejor} (validación)")
plt.xlabel("Minutos hasta el colegio")
plt.ylabel("Densidad")
plt.title("KDE con tres anchos de banda")
# Solo para dibujar: espacio arriba para la leyenda
plt.ylim(0, 0.19)
plt.legend(fontsize=11, ncol=2)
plt.show()

# Los mismos centros y el mismo estiramiento, con otra semilla
X_m2, _ = make_blobs(n_samples=30, centers=centros, cluster_std=1.15,
                     random_state=1)
nuevos_normales = X_m2 @ estirar
raros = np.array([[-6.0, 1.0], [1.0, 3.5], [-1.0, -2.5], [-1.0, 2.0]])
nuevos = np.vstack([nuevos_normales, raros])

umbral = np.percentile(gmm.score_samples(X_est), 1)
anomalia = gmm.score_samples(nuevos) < umbral
print(f"Umbral de log-densidad: {umbral:.2f}")
print("Marcados entre los 30 normales:", anomalia[:30].sum())
print("Marcados entre los 4 raros:    ", anomalia[30:].sum())

log_densidad = gmm.score_samples(rejilla).reshape(A.shape)
fig, eje = plt.subplots()
eje.contour(A, B, log_densidad, levels=[umbral], colors="C0",
            linestyles="--", linewidths=1.2)
eje.scatter(X_est[:, 0], X_est[:, 1], s=8, color="C2", alpha=0.5,
            label="Entrenamiento")
eje.scatter(nuevos[~anomalia, 0], nuevos[~anomalia, 1], s=36,
            color="C0", label="Nuevo: normal")
eje.scatter(nuevos[anomalia, 0], nuevos[anomalia, 1], s=110, marker="X",
            color="C1", edgecolors="white", label="Nuevo: anomalía")
eje.set_title("Densidad baja = anomalía")
# Solo para dibujar: espacio a la derecha para la leyenda
eje.set_xlim(-7, 7)
eje.legend(loc="upper right", fontsize=11)
plt.show()

busqueda_2d = GridSearchCV(KernelDensity(),
                           {"bandwidth": [0.03, 0.05, 0.1, 0.2, 0.4]},
                           cv=KFold(5, shuffle=True, random_state=0))
kde_lunas = busqueda_2d.fit(X_lunas).best_estimator_
gmm_lunas = GaussianMixture(n_components=2, random_state=0).fit(X_lunas)
print("Ancho de banda elegido:", kde_lunas.bandwidth)

X_lunas_nuevas, _ = make_moons(n_samples=300, noise=0.08, random_state=1)
for nombre, m in [("Mezcla de 2", gmm_lunas), ("KDE", kde_lunas)]:
    media = m.score_samples(X_lunas_nuevas).mean()
    print(f"{nombre:12s} log-densidad media en datos nuevos: {media:.2f}")

A2, B2 = np.meshgrid(np.linspace(-1.5, 2.5, 200),
                     np.linspace(-1, 1.5, 200))
rejilla2 = np.column_stack([A2.ravel(), B2.ravel()])
fig, ejes = plt.subplots(2, 1, figsize=(4.8, 4.8), sharex=True)
for eje, (nombre, m) in zip(ejes, [("Mezcla de 2 gaussianas", gmm_lunas),
                                   ("KDE", kde_lunas)]):
    dens = np.exp(m.score_samples(rejilla2)).reshape(A2.shape)
    eje.contour(A2, B2, dens, levels=6, colors="C2", linewidths=0.8)
    eje.scatter(X_lunas[:, 0], X_lunas[:, 1], s=6, color="C0")
    eje.set_title(nombre)
plt.show()

from sklearn.cluster import SpectralClustering

metodos = {
    "K-means": KMeans(n_clusters=2, n_init=10, random_state=0),
    "DBSCAN": DBSCAN(eps=0.2, min_samples=5),
    "Jerárquico simple": AgglomerativeClustering(n_clusters=2,
                                                 linkage="single"),
    "Espectral": SpectralClustering(n_clusters=2,
                                    affinity="nearest_neighbors",
                                    random_state=0),
    "Mezcla gaussiana": GaussianMixture(n_components=2, random_state=0),
}

fig, ejes = plt.subplots(3, 2, figsize=(4.8, 4.8),
                         sharex=True, sharey=True)
for eje, (nombre, metodo) in zip(ejes.flat, metodos.items()):
    grupos = metodo.fit_predict(X_lunas)
    acuerdo = adjusted_rand_score(y_lunas, grupos)
    print(f"{nombre:18s} acuerdo con las lunas: {acuerdo:.2f}")
    dibujar_grupos(eje, X_lunas, grupos, s=8)
    eje.set_title(nombre)
dibujar_grupos(ejes[2, 1], X_lunas, y_lunas, s=8)
ejes[2, 1].set_title("Las lunas reales")
plt.show()

from sklearn.datasets import load_sample_image

foto = load_sample_image("china.jpg")
# Una fila por píxel, con sus tres canales entre 0 y 1
pixeles = foto.reshape(-1, 3) / 255
print("Píxeles:", pixeles.shape)
print("Colores distintos en la foto:", len(np.unique(pixeles, axis=0)))

rng = np.random.default_rng(0)
muestra = pixeles[rng.choice(len(pixeles), size=2000, replace=False)]

# Solo para dibujar: las fotos son apaisadas, así que la figura es más baja
fig, ejes = plt.subplots(2, 2, figsize=(4.8, 3.8))
ejes[0, 0].imshow(foto)
ejes[0, 0].set_title("Original")
for eje, k in zip(ejes.flat[1:], [2, 8, 32]):
    colores = KMeans(n_clusters=k, n_init=4, random_state=0).fit(muestra)
    # Cada píxel se reemplaza por el centroide (el color) de su grupo
    grupo = colores.predict(pixeles)
    reducida = colores.cluster_centers_[grupo].reshape(foto.shape)
    eje.imshow(reducida)
    eje.set_title(f"{k} colores")
for eje in ejes.flat:
    eje.axis("off")
plt.show()