3.3 Agrupamiento: clasificación no supervisada - Machine Learning en las Geociencias

🖥️ Diapositivas — Sesión 13 (mié 28 oct)

El agrupamiento (clustering) es una forma de clasificación no supervisada: el algoritmo descubre la estructura de los datos a partir de las características (variables explicativas) únicamente, sin etiquetas. La meta es agrupar las observaciones en subgrupos coherentes.

Que un grupo sea coherente depende de la distancia entre los puntos de datos. La métrica de distancia cuantifica qué tan similares o disímiles son dos puntos de datos, y todo algoritmo de agrupamiento depende de una.

Este tutorial no cubre todos los métodos de agrupamiento posibles. Ningún método de agrupamiento funciona mejor en todos los escenarios; la elección correcta depende fuertemente de la estructura inherente de los datos. Un resumen simple, con ejemplos de juguete de estructuras de datos en 2D, está disponible en el paquete sklearn.

Esta lección se centra en conceptos fundamentales relevantes para las geociencias: 1) la definición de distancia, ilustrada con el algoritmo de agrupamiento más popular, 2) el agrupamiento k-means y 3) el agrupamiento aglomerativo.

1. Distancia¶

La distancia es la medida básica de disimilitud entre puntos de datos, y los resultados del agrupamiento dependen directamente de la métrica que usted elija. La distancia vuelve más adelante en el curso como el bloque constructivo de las funciones de pérdida y de costo para entrenar modelos de aprendizaje profundo.

Hay varias maneras de estimar y cuantificar la distancia entre puntos de datos. Algunas de las métricas de distancia más usadas son:

Distancia euclidiana: es la distancia en línea recta entre dos puntos de datos en un espacio multidimensional. Se usa a menudo cuando las características de los datos tienen unidades o escalas similares. $d(\mathbf{x},\mathbf{y}) = \sqrt{ \sum_{i=1}^N (x_i-y_i)^2 }$
Distancia de Manhattan: también conocida como distancia ‘L1’, mide la suma de las diferencias absolutas entre los elementos correspondientes de dos puntos de datos. Es adecuada cuando el movimiento a lo largo de los ejes está restringido, como en datos sobre mallas. $d(\mathbf{x},\mathbf{y}) = \sum_{i=1}^N |x_i-y_i|$
Distancia geodésica: la distancia geodésica mide el camino más corto entre dos puntos sobre la superficie de una esfera. Importa para la geografía terrestre, incluida la navegación GPS y las mediciones geodésicas.
Distancias basadas en correlación: en el análisis de datos geofísicos y geoespaciales, las métricas de distancia basadas en correlación, como la correlación de Pearson o la correlación de rangos de Spearman, se usan a menudo para evaluar relaciones entre variables.
Distancias basadas en covarianza: estas distancias, que toman en cuenta la covarianza espacial o los modelos de variograma, son frecuentes en la geoestadística y el análisis espacial.
Similitud coseno: esta métrica calcula el coseno del ángulo entre dos vectores de datos, y provee una medida de su similitud, en particular en espacios de alta dimensión. Se usa con frecuencia para datos de texto o de imágenes.

Scikit-learn contiene las métricas más usadas en el contexto del ML clásico en el paquete metrics.DistanceMetric. Más detalles en esta documentación de scikit-learn.

Relación con PCA El agrupamiento y PCA simplifican los datos mediante un número pequeño de resúmenes. Pero las diferencias son:

PCA busca reducir la dimensionalidad de los datos, encontrar una representación de baja dimensión que explique una buena fracción de la varianza de los datos,
el agrupamiento busca encontrar grupos homogéneos dentro de las observaciones.

De hecho, es común combinar ambos para datos complejos y de alta dimensión: 1) PCA, 2) agrupamiento sobre las componentes principales.

Hay dos métodos principales de agrupamiento: el agrupamiento k-means y el agrupamiento jerárquico.

La caja de herramientas scikit-learn tiene una colección de algoritmos de agrupamiento y una documentación detallada con tutoriales.

2. Preparación del tutorial¶

Importar los paquetes de Python útiles.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pooch
from sklearn import preprocessing
from sklearn.decomposition import PCA

# one generator for the whole notebook, the idiom from 2.6
rng = np.random.default_rng(42)

Datos¶

Las ediciones anteriores de este capítulo usaban una tabla de citometría de flujo de SeaFlow de 146 MB; la retiramos en 2026 en favor de conjuntos de datos más ligeros.

Nuestro primer conjunto de datos de trabajo es el registro del géiser Old Faithful, en Yellowstone. Cada fila empareja la duración de una erupción (current, en minutos) con la duración de la erupción siguiente (next, en minutos). Lo descargamos con pooch, que guarda el archivo en una caché local.

path = pooch.retrieve(
    url="https://raw.githubusercontent.com/UW-MLGEO/MLGeo-dataset/main/data/faithful.csv",
    known_hash=None,
    fname="faithful.csv",
    path=pooch.os_cache("mlgeo"),
)
faithful = pd.read_csv(path)
faithful.head()

plt.plot(faithful.current, faithful.next, 'ko')
plt.xlabel('Eruption duration in minutes')
plt.ylabel('Duration of next eruption in minutes')
plt.grid(True)
plt.title('Old Faithful')

K-means depende de la distancia euclidiana, así que las características con rangos numéricos mayores dominan el resultado. Aquí ambas columnas están en minutos con rangos similares, pero eso es una coincidencia de este conjunto de datos. Estandarizar cada característica (media cero, varianza uno) es un hábito que conviene conservar para cualquier método basado en distancias — en el ejercicio sísmico del final de este cuaderno, las características abarcan órdenes de magnitud y el escalado no es opcional. Estandarizamos ahora y guardamos el resultado en data_faithful, un arreglo de numpy de 2 columnas.

scaler = preprocessing.StandardScaler()
data_faithful = scaler.fit_transform(faithful[["current", "next"]])
print(data_faithful.shape)

(271, 2)

3. K-means¶

K-means (k-medias) es un método de agrupamiento no supervisado. La idea principal es separar los datos en K conglomerados (clusters) distintos. Para los grupos que devuelve un algoritmo, «conglomerado», «grupo» y «clúster» circulan por igual según la comunidad; este cuaderno dice «conglomerado» de forma consistente, y «clúster» queda reservado para la infraestructura de cómputo (capítulo 5). Tenemos entonces dos problemas que resolver. Primero, hay que encontrar los k centroides de los k conglomerados. Luego, hay que asignar cada punto de datos al conglomerado cuyo centroide le queda más cerca.

La meta es particionar $n$ puntos de datos en $k$ conglomerados. Cada observación se etiqueta con el conglomerado de media más cercana.

K-means es iterativo:

suponga valores iniciales para la media de cada uno de los $k$ conglomerados
calcule la distancia de cada observación a cada una de las $k$ medias
etiquete cada observación como perteneciente a la media más cercana
encuentre el centro de masa (la media) de cada grupo de puntos etiquetados. Estas son las nuevas medias para el paso 1.

En lo que sigue, denotamos $n$ el número de puntos de datos y $p$ el número de características de cada punto de datos.

K-means solo funciona con la métrica de distancia euclidiana. De hecho, su concepto central es usar la distancia euclidiana para medir y minimizar la inercia o variación intraconglomerado.

n, p = data_faithful.shape
print('We have {:d} data points, and each one has {:d} features'.format(n, p))

We have 271 data points, and each one has 2 features

Definamos una función para inicializar los centroides de los conglomerados. Elegimos puntos aleatorios dentro del rango de valores que toman los datos.

def init_centers(data, k):
    """
    """
    # Initialize centroids
    centers = np.zeros((k, np.shape(data)[1]))
    # Loop on k centers
    for i in range(0, k):
        # Generate p random values between 0 and 1
        dist = rng.uniform(size=np.shape(data)[1])
        # Use the random values to generate a point within the range of values taken by the data
        centers[i, :] = np.min(data, axis=0) + (np.max(data, axis=0) - np.min(data, axis=0)) * dist
    return centers

Para poder asignar cada punto de datos al centroide más cercano, necesitamos definir la distancia entre dos puntos de datos. La distancia más común es la distancia euclidiana:

$d(x,y) = \sqrt{\sum_{i = 1}^p (x_i - y_i)^2}$

donde $x$ y $y$ son dos puntos de observación con $p$ variables.

Definimos entonces una función para calcular la distancia entre cada punto de datos y cada centroide.

def compute_distance(data, centers, k):
    """
    """
    # Initialize distance
    distance = np.zeros((np.shape(data)[0], k))
    # Loop on n data points
    for i in range(0, np.shape(data)[0]):
        # Loop on k centroids
        for j in range(0, k):
            # Compute distance
            distance[i, j] = np.sqrt(np.sum(np.square(data[i, :] - centers[j, :])))
    return distance

Ahora definimos una función para asignar cada punto de datos al conglomerado cuyo centroide le queda más cerca. También definimos una función objetivo que se minimizará hasta alcanzar la convergencia.

Nuestro objetivo es minimizar la suma de los cuadrados de la distancia entre cada punto y el centroide más cercano:

$obj = \sum_{j = 1}^k \sum_{i = 1}^{N_j} d(x^{(i)} , x^{(j)}) ^2$

donde $x^{(i)}$ es el $i$ -ésimo punto del conglomerado $j$ , $x^{(j)}$ es el centroide del conglomerado $j$ , y $N_j$ es el número de puntos del conglomerado $j$ .

def compute_objective(distance, clusters):
    """
    """
    # Initialize objective
    objective = 0.0
    # Loop on n data points
    for i in range(0, np.shape(distance)[0]):
        # Add distance to the closest centroid
        objective = objective + distance[i, int(clusters[i])] ** 2.0
    return objective

def compute_clusters(distance):
    """
    """
    # Initialize clusters
    clusters = np.zeros(np.shape(distance)[0])
    # Loop on n data points
    for i in range(0, np.shape(distance)[0]):
        # Find closest centroid
        best = np.argmin(distance[i, :])
        # Assign data point to corresponding cluster
        clusters[i] = best
    return clusters

Después de asignar todos los puntos a un conglomerado, calcule la nueva ubicación del centroide. Es simplemente el valor de la media de todos los puntos asignados a ese conglomerado:

Para $1 \leq j \leq k$ , $x_p^{(j)} = \frac{1}{N_j} \sum_{i = 1}^{N_j} x_p^{(i)}$

def compute_centers(data, clusters, k):
    """
    """
    # Initialize centroids
    centers = np.zeros((k, np.shape(data)[1]))
    # Loop on clusters
    for i in range(0, k):
        # Select all data points in this cluster
        subdata = data[clusters == i, :]
        # If no data point in this cluster, generate randomly a new centroid
        if (np.shape(subdata)[0] == 0):
            centers[i, :] = init_centers(data, 1)
        else:
            # Compute the mean location of all data points in this cluster
            centers[i, :] = np.mean(subdata, axis=0)
    return centers

Ya podemos programar el algoritmo k-means ensamblando todas estas funciones. Detenemos el cálculo cuando la función objetivo deja de disminuir.

def my_kmeans(data, k):
    """
    """
    # Initialize centroids
    centers = init_centers(data, k)
    # Initialize objective function to square of the maximum distance between two data points times number of data points
    objective_old = np.shape(data)[0] * np.sum(np.square(np.max(data, axis=0) - np.min(data, axis=0)))
    # Initialize clusters
    clusters_old = np.zeros(np.shape(data)[0])
    # Start loop until convergence
    stop_alg = False
    while stop_alg == False:
        # Compute distance between data points and centroids
        distance = compute_distance(data, centers, k)
        # Get new clusters
        clusters_new = compute_clusters(distance)
        # get new value of objective function
        objective_new = compute_objective(distance, clusters_new)
        # If objective function stops decreasing, end loop
        if objective_new >= objective_old:
            return (clusters_old, objective_old, centers)
        else:
            # Update the locations of the centroids
            centers = compute_centers(data, clusters_new, k)
            objective_old = objective_new
            clusters_old = clusters_new

Los datos de Old Faithful muestran dos regímenes — erupciones cortas seguidas de erupciones cortas, y erupciones largas seguidas de erupciones largas —, así que fijamos k = 2 y corremos nuestro k-means hecho desde cero.

k = 2
(clusters, objective, centers) = my_kmeans(data_faithful, k)

plt.scatter(data_faithful[:, 0], data_faithful[:, 1], c=clusters)
plt.scatter(centers[:, 0], centers[:, 1], marker='o', s=300, c='black')
plt.xlabel('Eruption duration (standardized)')
plt.ylabel('Next eruption duration (standardized)')
plt.grid(True)
plt.title('From-scratch k-means with k = 2')

K-means con Scikit-learn¶

Ahora usamos la implementación de scikit-learn sobre los mismos datos escalados.

Dos valores por defecto que vale la pena conocer: init='k-means++' es el esquema de inicialización por defecto, y desde las versiones recientes de scikit-learn n_init toma por defecto 'auto', que corre varios reinicios de k-means++ y conserva el mejor. Ya no hace falta fijar ninguno de los dos parámetros a mano.

from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=2, random_state=0).fit(data_faithful)

plt.scatter(data_faithful[:, 0], data_faithful[:, 1], c=kmeans.labels_)
plt.scatter(kmeans.cluster_centers_[:, 0], kmeans.cluster_centers_[:, 1], marker='o', s=300, c='black')
plt.xlabel('Eruption duration (standardized)')
plt.ylabel('Next eruption duration (standardized)')
plt.grid(True)
plt.title('Scikit-learn KMeans with k = 2')

Consejos prácticos para k-means¶

¿Cómo evaluar el éxito del agrupamiento? ¿Qué tan bien separados están los conglomerados? Existen muchas maneras de evaluar la calidad de los conglomerados. Scikit-learn ha resumido y empaquetado estas herramientas en el módulo sklearn.metrics. Vea la documentación aquí. En general, los scores altos son mejores. Aquí un resumen:
- Si los datos tienen etiquetas verdaderas (ground truth; por ejemplo, los tipos de fuente sísmica del ejercicio de abajo), podemos usar varias métricas:
  - homogeneidad (cada conglomerado contiene solo miembros de una clase dada, metrics.homogeneity_score(clusterID,true_label)), completitud (todos los miembros de una clase dada se asignan al mismo conglomerado, metrics.completeness_score), la medida V (2 x homogeneidad x completitud / (homogeneidad+completitud), metrics.v_measure_score) y las tres juntas metrics.homogeneity_completeness_v_measure(clusterID,true_label).
  - el índice de Fowlkes-Mallows, FMI, que usa TP (verdaderos positivos), FP (falsos positivos) y FN (falsos negativos). Vale 0.0 para una asignación aleatoria de conglomerados y 1.0 para asignaciones de etiquetas perfectas.
- Si los datos no tienen etiquetas verdaderas, puede usar:
  - el coeficiente de silueta con el módulo metrics.silhouette_score; un score o coeficiente alto es mejor. Cuantifica qué tan compactos son los datos dentro de los conglomerados y qué tan bien separados están los conglomerados. Más detalles sobre la implementación y visualización del score de silueta aquí.
El número de conglomerados $k$ es un parámetro ajustable. Para encontrar el número óptimo de conglomerados, discutimos abajo algunas estrategias.
La optimización de k-means puede tener mínimos locales. El resultado, por lo tanto, puede diferir según la inicialización de los conglomerados. Se recomienda:
- repetir la inicialización aleatoria, repetir k-means y usar el mejor conjunto de conglomerados (el que tenga el menor error final)
- elegir el esquema de inicialización K-means++ con el parámetro de sklearn init='k-means++'. El algoritmo selecciona los centroides iniciales mediante un muestreo basado en una distribución de probabilidad empírica de la contribución de los puntos a la inercia total.
La varianza de algunas características puede afectar los resultados. Puede ser difícil encontrar conglomerados si algunas características de los datos (ejes) son mucho más grandes que otras. Los datos pueden requerir preprocesamiento (centrado y escalado) o preacondicionamiento (por ejemplo, PCA).

Análisis de silueta

El coeficiente de silueta y el score de silueta (silhouette score) son métricas para evaluar la calidad de un agrupamiento en el aprendizaje no supervisado.

Coeficiente de silueta: El coeficiente de silueta de una muestra de datos mide qué tan similar es a su propio conglomerado (cohesión) en comparación con los otros conglomerados (separación). Se calcula para cada muestra de datos y va de -1 a 1. Un coeficiente de silueta alto indica que el objeto está bien emparejado con su propio conglomerado y mal emparejado con los conglomerados vecinos. A la inversa, un coeficiente de silueta bajo sugiere que el objeto puede estar en el conglomerado equivocado.
La fórmula del coeficiente de silueta (s) para un punto de datos individual es:
$s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}}$ donde:
- $a(i)$ es la distancia promedio del $i$ -ésimo punto de datos a los demás puntos de su mismo conglomerado.
- $b(i)$ es la menor distancia promedio del $i$ -ésimo punto de datos a los puntos de un conglomerado distinto, minimizada sobre los conglomerados.
Score de silueta: El score de silueta es el promedio del coeficiente de silueta sobre todas las muestras de datos del conjunto. Provee una medida global de qué tan bien separados están los conglomerados. El score de silueta también va de -1 a 1: un score alto indica un buen agrupamiento y un score bajo sugiere conglomerados superpuestos o mal clasificados.
La fórmula del score de silueta es:

$S = \frac{\sum_{i=1}^{N} s(i)}{N}$ donde:

$N$ es el número de puntos de datos del conjunto.

Interpretación:

Un coeficiente de silueta cercano a 1 indica que el punto de datos está bien emparejado con su propio conglomerado y mal emparejado con los conglomerados vecinos, lo que señala un conglomerado robusto y bien diferenciado.
Un coeficiente de silueta cercano a -1 sugiere que el punto de datos posiblemente está mal clasificado, pues empareja mejor con un conglomerado vecino.
Un coeficiente de silueta alrededor de 0 indica conglomerados superpuestos.

Para el score de silueta:

Un score cercano a 1 implica conglomerados bien definidos.
Un score alrededor de 0 sugiere conglomerados superpuestos.
Un score negativo indica que la mayoría de los puntos de datos puede estar asignada a los conglomerados equivocados.

En resumen, el análisis de silueta ayuda a seleccionar el número de conglomerados y a evaluar la calidad general de un agrupamiento.

# example of the silhouette score for the Old Faithful data
from sklearn.cluster import KMeans
from sklearn import metrics
from sklearn.metrics import silhouette_score, silhouette_samples

ncluster = 2

kmeans_model = KMeans(n_clusters=ncluster, random_state=1).fit(data_faithful)
labels = kmeans_model.labels_
sc = silhouette_score(data_faithful, labels, metric='euclidean')
print("Silhouette score for k = 2:", sc)

Silhouette score for k = 2: 0.5611599583255005

import matplotlib.cm as cm
fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 7)
# The 1st subplot is the silhouette plot
# The silhouette coefficient can range from -1, 1 but in this example all
# lie within [-0.1, 1]
ax1.set_xlim([-0.1, 1])
# The (n_clusters+1)*10 is for inserting blank space between silhouette
# plots of individual clusters, to demarcate them clearly.
ax1.set_ylim([0, len(data_faithful) + (ncluster + 1) * 10])

# Initialize the clusterer with n_clusters value and a random generator
# seed of 10 for reproducibility.
clusterer = KMeans(n_clusters=ncluster, random_state=10)
cluster_labels = clusterer.fit_predict(data_faithful)

# The silhouette_score gives the average value for all the samples.
# This gives a perspective into the density and separation of the formed
# clusters
silhouette_avg = silhouette_score(data_faithful, cluster_labels)
print(
    "For n_clusters =",
    ncluster,
    "The average silhouette_score is :",
    silhouette_avg,
)

# Compute the silhouette scores for each sample
sample_silhouette_values = silhouette_samples(data_faithful, cluster_labels)

y_lower = 10
for i in range(ncluster):
    # Aggregate the silhouette scores for samples belonging to
    # cluster i, and sort them
    ith_cluster_silhouette_values = sample_silhouette_values[cluster_labels == i]

    ith_cluster_silhouette_values.sort()

    size_cluster_i = ith_cluster_silhouette_values.shape[0]
    y_upper = y_lower + size_cluster_i

    color = cm.nipy_spectral(float(i) / ncluster)
    ax1.fill_betweenx(
        np.arange(y_lower, y_upper),
        0,
        ith_cluster_silhouette_values,
        facecolor=color,
        edgecolor=color,
        alpha=0.7,
    )

    # Label the silhouette plots with their cluster numbers at the middle
    ax1.text(-0.05, y_lower + 0.5 * size_cluster_i, str(i))

    # Compute the new y_lower for next plot
    y_lower = y_upper + 10  # 10 for the 0 samples

ax1.set_title("The silhouette plot for the various clusters.")
ax1.set_xlabel("The silhouette coefficient values")
ax1.set_ylabel("Cluster label")

# The vertical line for average silhouette score of all the values
ax1.axvline(x=silhouette_avg, color="red", linestyle="--")

ax1.set_yticks([])  # Clear the yaxis labels / ticks
ax1.set_xticks([-0.1, 0, 0.2, 0.4, 0.6, 0.8, 1])

# 2nd Plot showing the actual clusters formed
colors = cm.nipy_spectral(cluster_labels.astype(float) / ncluster)
ax2.scatter(
    data_faithful[:, 0], data_faithful[:, 1], marker="o", s=30, lw=0, alpha=0.7, c=colors, edgecolor="k"
)

# Labeling the clusters
centers = clusterer.cluster_centers_
# Draw white circles at cluster centers
ax2.scatter(
    centers[:, 0],
    centers[:, 1],
    marker="o",
    c="white",
    alpha=1,
    s=200,
    edgecolor="k",
)

for i, c in enumerate(centers):
    ax2.scatter(c[0], c[1], marker="$%d$" % i, alpha=1, s=200, edgecolor=None)

ax2.set_title("The visualization of the clustered data.")
ax2.set_xlabel("Feature space for the 1st feature")
ax2.set_ylabel("Feature space for the 2nd feature")

plt.suptitle(
    "Silhouette analysis for KMeans clustering on sample data with n_clusters = %d"
    % ncluster,
    fontsize=14,
    fontweight="bold",
)

For n_clusters = 2 The average silhouette_score is : 0.5611599583255005

Elección del número de conglomerados: el método del codo¶

El método del codo está diseñado para encontrar el número óptimo de conglomerados. Consiste en ejecutar el algoritmo de agrupamiento con un número creciente de conglomerados $k$ , midiendo la distancia promedio entre los puntos de datos y los centroides. Hay dos métricas típicas en el método del codo:

Distorsión: se calcula como el promedio de las distancias al cuadrado a los centros de los conglomerados respectivos. Típicamente se usa la métrica de distancia euclidiana.
Inercia: es la suma de las distancias al cuadrado de las muestras a su centro de conglomerado más cercano.

Para cada valor de $k$ , calculamos la media del cuadrado de la distancia entre los puntos de datos y el centroide del conglomerado al que pertenecen. Luego graficamos ese valor en función de $k$ . Con suerte, disminuye y luego alcanza una meseta. El número óptimo de conglomerados es el valor en el quiebre — el «codo» de la curva.

Ilustramos el método sobre los datos estandarizados de Old Faithful.

def compute_elbow(data, clusters, centers, k):
    """
    """
    E = 0
    for i in range(0, k):
        distance = compute_distance(data[clusters == i, :], centers[i, :].reshape(1, -1), 1)
        E = E + np.mean(np.square(distance))
    return E

Calcule el valor de E para k entre 1 y 8 y grafíquelo. Busque el valor de $k$ donde la curva se dobla.

E = np.zeros(8)
for k in range(1, 9):
    (clusters, objective, centers) = my_kmeans(data_faithful, k)
    E[k - 1] = compute_elbow(data_faithful, clusters, centers, k)

plt.figure(figsize=(4, 4))
plt.plot(np.arange(1, 9), E)
plt.xlabel('Number of clusters')
plt.ylabel('Elbow criterion')
plt.grid(True)

El método del codo no siempre funciona bien. Por ejemplo, vea qué pasa cuando los puntos se acercan entre sí. Como data_faithful está estandarizado y centrado en cero, encogemos los datos hacia el origen.

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 7)
alpha = 0.5
origin = np.array([0.0, 0.0])
data_shrink = origin + alpha * np.sign(data_faithful - origin) * np.power(np.abs(data_faithful - origin), 2.0)
ax1.plot(data_shrink[:, 0], data_shrink[:, 1], 'ko')
ax1.set_title('Shrunken data')

E = np.zeros(8)
for k in range(1, 9):
    (clusters, objective, centers) = my_kmeans(data_shrink, k)
    E[k - 1] = compute_elbow(data_shrink, clusters, centers, k)
ax2.plot(np.arange(1, 9), E)
ax2.set_xlabel('Number of clusters')
ax2.set_ylabel('Elbow criterion')
ax2.grid(True)

Veamos qué pasa cuando disminuimos el número de datos.

fig, (ax1, ax2) = plt.subplots(1, 2)
fig.set_size_inches(18, 7)
alpha = 0.2
indices = rng.uniform(size=np.shape(data_faithful)[0])
subdata = data_faithful[indices < alpha, :]
ax1.plot(subdata[:, 0], subdata[:, 1], 'ko')
ax1.set_title('Subsampled data')
E = np.zeros(8)
for k in range(1, 9):
    (clusters, objective, centers) = my_kmeans(subdata, k)
    E[k - 1] = compute_elbow(subdata, clusters, centers, k)
ax2.plot(np.arange(1, 9), E)
ax2.set_xlabel('Number of clusters')
ax2.set_ylabel('Elbow criterion')
ax2.grid(True)

Ambas pruebas de estrés muestran los límites del criterio del codo: cuando los conglomerados se acercan entre sí, o cuando los datos se vuelven escasos, la curva pierde su quiebre claro y la elección de $k$ se vuelve ambigua. Use el método del codo como guía, no como regla, y contrástelo con el score de silueta.

Repetir k-means¶

El resultado es muy sensible a la ubicación de los centroides iniciales. Repita el agrupamiento N veces y elija el agrupamiento con la mejor función objetivo.

def repeat_kmeans(data, k, N):
    """
    """
    # Initialization
    objective = np.zeros(N)
    clusters = np.zeros((N, np.shape(data)[0]))
    centers = np.zeros((N, k, np.shape(data)[1]))
    # Run K-means N times
    for i in range(0, N):
        result = my_kmeans(data, k)
        clusters[i, :] = result[0]
        objective[i] = result[1]
        centers[i, :, :] = result[2]
    # Choose the clustering with the best value of the objective function
    best = np.argmin(objective)
    return (clusters[best, :], objective[best], centers[best, :, :])

Repita k-means 20 veces con k = 2 y conserve la mejor corrida.

N = 20
k = 2
(clusters, objective, centers) = repeat_kmeans(data_faithful, k, N)

plt.figure(figsize=(6, 6))
plt.scatter(data_faithful[:, 0], data_faithful[:, 1], c=clusters)
plt.scatter(centers[:, 0], centers[:, 1], marker='o', s=300, c='black')
plt.xlabel('Eruption duration (standardized)')
plt.ylabel('Next eruption duration (standardized)')
plt.grid(True)
plt.title('Best of {:d} k-means runs'.format(N))

K-means puede ser un proceso lento de calcular y hay vías para acelerarlo. Una solución es usar k-means por mini-lotes (mini-batch), que toma un subconjunto de los datos en cada iteración para construir los centroides.

4. Agrupamiento jerárquico¶

En k-means, usamos la distancia euclidiana y prescribimos el número de conglomerados K.

En el agrupamiento jerárquico, elegimos distintas métricas de distancia, visualizamos la estructura de los datos y luego decidimos el número de conglomerados. Hay dos enfoques para construir la jerarquía de conglomerados:

Aglomerativo: cada punto comienza en su propio conglomerado. Los datos se fusionan por pares mientras se crea una jerarquía de conglomerados.
Divisivo: al inicio, todos los datos están en 1 conglomerado. Los datos se dividen recursivamente en conglomerados cada vez más pequeños.

Hay varios tipos de enlaces (linkages). sklearn tiene documentación detallada, sobre todo para el enfoque aglomerativo. Los distintos métodos de enlace son:

Ward minimiza la suma de las diferencias al cuadrado dentro de todos los conglomerados. Es un enfoque de minimización de varianza y, en ese sentido, es similar a la función objetivo de k-means, pero abordado con un enfoque jerárquico aglomerativo.
El enlace máximo o completo minimiza la distancia máxima entre observaciones de pares de conglomerados.
El enlace promedio minimiza el promedio de las distancias entre todas las observaciones de pares de conglomerados.
El enlace simple minimiza la distancia entre las observaciones más cercanas de pares de conglomerados.

Primero importamos los paquetes pertinentes.

from matplotlib import rcParams
from scipy.cluster import hierarchy
from scipy.spatial.distance import pdist

rcParams.update({'font.size': 18})
plt.rcParams['figure.figsize'] = [12, 12]

Aquí creamos un conjunto de datos ficticio con 2 conglomerados que se entremezclan en unos pocos puntos de datos.

# Training and testing set sizes
n = 100 # Train

# Random ellipse 1 centered at (0,0)
x = rng.standard_normal(n)
y = 0.5*rng.standard_normal(n)

# Random ellipse 2 centered at (1,-2)
x2 = rng.standard_normal(n) + 1
y2 = 0.2*rng.standard_normal(n) - 2

# Rotate ellipse 2 by theta
theta = np.pi/4
A = np.zeros((2,2))
A[0,0] = np.cos(theta)
A[0,1] = -np.sin(theta)
A[1,0] = np.sin(theta)
A[1,1] = np.cos(theta)

x3 = A[0,0]*x2 + A[0,1]*y2
y3 = A[1,0]*x2 + A[1,1]*y2

plt.figure(1,figsize=(5,5))
plt.plot(x[:],y[:],'ro')
plt.plot(x3[:],y3[:],'bo')
plt.grid(True)
plt.show()

Combinar estos dos conjuntos de datos en uno.

X1 = np.column_stack((x3[:],y3[:]))
X2 = np.column_stack((x[:],y[:]))
X = np.concatenate((X1,X2))

plt.figure(1,figsize=(5,5))
plt.plot(X[:,0],X[:,1],'ro')
plt.grid(True)
plt.show()

Primero exploramos los dendrogramas.

## Dendrograms
Y = pdist(X,metric='euclidean')
Z = hierarchy.linkage(Y,method='average')
thresh = 0.5*np.max(Z[:,2])

plt.figure(1,figsize=(5,5))
dn = hierarchy.dendrogram(Z,p=100,color_threshold=thresh)
plt.xlabel('Data Sample Index')
plt.ylabel('Distance')
plt.title('Dendrogram with average linkage, high threshold')
plt.show()

thresh = 0.25*np.max(Z[:,2])
plt.figure()
dn = hierarchy.dendrogram(Z,p=100,color_threshold=thresh)
plt.xlabel('Data Sample Index')
plt.ylabel('Distance')
plt.title('Dendrogram with average linkage, low threshold')
plt.show()

Ya exploramos la estructura de los datos y construimos algo de intuición sobre cuántos conglomerados hay y cómo se distribuyen las muestras de datos entre ellos.

A continuación, elegimos un umbral de distancia y asignamos a cada punto de datos un identificador de conglomerado.

En Scikit-learn, el algoritmo completo está incorporado en la función AgglomerativeClustering documentación de sklearn aquí. La función necesita de todos modos un umbral de distancia o un número de conglomerados para realizar el agrupamiento y asignar un identificador de conglomerado a cada muestra de datos.

from sklearn.cluster import AgglomerativeClustering
# Let's first find a reasonable distance threshold by precalculating the linkage matrix
Z = hierarchy.linkage(X, "average")
thresh = 0.85*np.max(Z[:,2])    # choose a threshold distance
print(thresh)
# design model
model = AgglomerativeClustering(distance_threshold=thresh,linkage="average", n_clusters=None)
# fit model and predict clusters on the data samples
clusterID = model.fit_predict(X)

3.349805679694126

plt.figure(figsize=(6,6))
plt.scatter(X[:,0],X[:,1],c=clusterID)

5. PCA antes del agrupamiento¶

Generemos datos sintéticos.

centers = np.array([[2, 2], [2, 8], [4, 3]])
radius = [0.1, 1]
synthetics = np.empty([0, 2])
for i in range(0, 3):
    X = centers[i, 0] + radius[0] * rng.standard_normal(100)
    Y = centers[i, 1] + radius[1] * rng.standard_normal(100)
    U = (X + Y) * np.sqrt(2) / 2
    V = (X - Y) * np.sqrt(2) / 2
    synthetics = np.concatenate([synthetics, np.vstack((U, V)).T])
plt.figure(figsize=(6,6))
plt.plot(synthetics[:,0], synthetics[:,1], 'ko')
plt.xlim(1, 9)
plt.ylim(-6, 3)

(-6.0, 3.0)

Hagamos ahora un agrupamiento k-means con 3 conglomerados.

(clusters, objective, centers) = my_kmeans(synthetics, 3)
plt.figure(figsize=(6,6))
plt.scatter(synthetics[:,0], synthetics[:,1], c=clusters)

¿Qué pasa si aplicamos PCA + normalización antes del agrupamiento?

pca = PCA(n_components=2)
synthetics_pca = pca.fit_transform(synthetics)

scaler = preprocessing.StandardScaler().fit(synthetics_pca)
synthetics_scaled = scaler.transform(synthetics_pca)

(clusters, objective, centers) = my_kmeans(synthetics_scaled, 3)
plt.figure(figsize=(6,6))
plt.scatter(synthetics[:,0], synthetics[:,1], c=clusters)

Ejercicio: agrupar fuentes sísmicas volcánicas y tectónicas¶

Los volcanes vigilados de cerca — el Popocatépetl, que monitorea el CENAPRED; los volcanes glaciados de los Andes, que monitorean los observatorios de Chile, Ecuador y Colombia; o el Mt Hood, un volcán cubierto de hielo en las Cascadas — producen una mezcla de fuentes sísmicas: sismos volcánicos y tectónicos, eventos superficiales como caídas de rocas y avalanchas, y ruido de fondo. Los analistas etiquetan estos eventos a mano. ¿Puede el agrupamiento encontrar los tipos de fuente a partir de las características de la forma de onda únicamente?

Usamos un conjunto de datos curado de eventos sísmicos del Noroeste del Pacífico estadounidense — sismo, explosión, evento superficial y ruido, 1000 de cada uno —, descritos por características físicas de la forma de onda (forma espectral, estadísticas de la envolvente, curtosis, energías por banda). El conjunto de datos está archivado en Zenodo: DOI 10.5281/zenodo.14025693. Regresa en el cuaderno 3.5 para la clasificación supervisada.

Paso 1: cargar los datos. El cargador de abajo descarga y guarda en caché los cuatro archivos de clase, los concatena y elimina la única columna de características con valores faltantes.

import pooch

SEISMIC_FILES = {
    "1000_earthquakes_physical_features.csv": "md5:28129c8dd1b3e14f655d489577b841b5",
    "1000_explosion_physical_features.csv": "md5:af1342d32e163e961e043364136359b0",
    "1000_noise_physical_features.csv": "md5:16cdb992fed6cf6273d5624f5df905da",
    "1000_surface_physical_features.csv": "md5:9a2c2643030cf058704d68e130654e9d",
}

frames = []
for fname, checksum in SEISMIC_FILES.items():
    path = pooch.retrieve(
        url=f"https://zenodo.org/api/records/14025693/files/{fname}/content",
        known_hash=checksum,
        fname=fname,
        path=pooch.os_cache("mlgeo"),
    )
    frames.append(pd.read_csv(path, index_col=0))
seismic = pd.concat(frames, ignore_index=True)
seismic = seismic.dropna(axis=1)  # drops the one feature column with missing values

X_seis = seismic.drop(columns=["source", "serial_no"])
y_seis = seismic["source"]
print(X_seis.shape)
y_seis.value_counts()

(4000, 61)

source
earthquake       1000
explosion        1000
noise            1000
surface event    1000
Name: count, dtype: int64

Paso 2: estandarizar las 61 características. Las características abarcan órdenes de magnitud, así que el escalado es obligatorio antes de cualquier método basado en distancias.

scaler = preprocessing.StandardScaler()
X_seis_scaled = scaler.fit_transform(X_seis)

Paso 3: PCA. Grafique la varianza explicada acumulada y conserve suficientes componentes principales para explicar alrededor del 80 % de la varianza.

pca = PCA().fit(X_seis_scaled)
cumvar = np.cumsum(pca.explained_variance_ratio_)

plt.figure(figsize=(6, 4))
plt.plot(np.arange(1, len(cumvar) + 1), cumvar, marker='.')
plt.axhline(0.80, color='red', linestyle='--', label='80% of variance')
plt.xlabel('Number of principal components')
plt.ylabel('Cumulative explained variance')
plt.legend()
plt.grid(True)

n_pcs = int(np.searchsorted(cumvar, 0.80) + 1)
print(f"{n_pcs} components explain 80% of the variance")
X_seis_pca = pca.transform(X_seis_scaled)[:, :n_pcs]

15 components explain 80% of the variance

Paso 4: k-means con 4 conglomerados sobre las componentes principales. Sabemos que hay cuatro tipos de fuente, así que fijamos k = 4.

kmeans_seis = KMeans(n_clusters=4, random_state=42).fit(X_seis_pca)
cluster_labels = kmeans_seis.labels_

Paso 5: comparar los conglomerados con las etiquetas verdaderas. Este conjunto de datos tiene etiquetas verdaderas (ground truth), así que podemos calificar el agrupamiento con tres métricas:

La homogeneidad vale 1 cuando cada conglomerado contiene miembros de una sola clase.
La completitud vale 1 cuando todos los miembros de una clase caen en el mismo conglomerado.
La medida V es la media armónica de las dos. Las tres van de 0 (asignación aleatoria) a 1 (correspondencia perfecta).

h, c, v = metrics.homogeneity_completeness_v_measure(y_seis, cluster_labels)
print(f"Homogeneity:  {h:.3f}")
print(f"Completeness: {c:.3f}")
print(f"V-measure:    {v:.3f}")

Homogeneity:  0.363
Completeness: 0.393
V-measure:    0.377

Paso 6: cruzar en una tabla los conglomerados contra los tipos de fuente.

ct = pd.crosstab(cluster_labels, y_seis, rownames=['cluster'], colnames=['source'])
ct

fig, ax = plt.subplots(figsize=(7, 5))
im = ax.imshow(ct.values, cmap='Blues')
ax.set_xticks(range(len(ct.columns)))
ax.set_xticklabels(ct.columns, rotation=45, ha='right')
ax.set_yticks(range(len(ct.index)))
ax.set_yticklabels([f'cluster {i}' for i in ct.index])
for i in range(ct.shape[0]):
    for j in range(ct.shape[1]):
        ax.text(j, i, ct.values[i, j], ha='center', va='center',
                color='white' if ct.values[i, j] > ct.values.max()/2 else 'black')
fig.colorbar(im, ax=ax, label='Number of events')
ax.set_title('K-means clusters vs. true source types')

Sin ver jamás una etiqueta, k-means recupera buena parte de la estructura: un conglomerado está dominado por los sismos y otro por el ruido. Mire en su tabla cruzada el par de fuentes que más se mezcla — aquí, las explosiones y los eventos superficiales comparten en gran medida un conglomerado. Eso es físicamente plausible: ambas son fuentes someras, cercanas a la superficie (explosiones de cantera, caídas de rocas, avalanchas), así que sus formas de onda llevan poca de la energía profunda de alta frecuencia que separa a los sismos, y sus características de envolvente y espectrales se superponen. Los clasificadores supervisados del cuaderno 3.5 lo hacen mucho mejor sobre las mismas características — ese es precisamente el aporte de las etiquetas.

References¶

Kharita, A. (2024). Physical features for small sample of data (1000 events per class). Zenodo. 10.5281/ZENODO.14025693