Optimizar el número de clústeres con gap statistics

El mayor problema con el que nos podemos encontrar a la hora de usar el algoritmo de k-means es conocer el número de clústeres en los que se divide el conjunto de datos. Un hiperparámetro que en Scikit-learn debe ser indicado al construir el objeto. Por eso existen múltiples métodos para seleccionar este valor como los métodos del codo (elbow method), silueta (silhouette) o índice de Calinski-Harabasz. Otro método bastante popular es el conocido como estadísticas del gap (gap statistics). Veamos en qué consiste y cómo se puede implementar en Python.

El método gap statistics

La forma que usa el método de gap statistics para estimar el número óptimo de clústeres es compara la dispersión dentro del conjunto de datos original con la dispersión que se puede esperar en conjuntos de datos generados aleatoriamente. Así, si se busca el punto en el que se maximiza la diferencias entre ambos valores se obtiene se puede estimar la cantidad óptima de clústeres.

El método gap statistics para obtener el número óptimo de clústeres se puede resumir en los siguientes pasos:

1. Seleccionar el número máximo de clústeres (k_{max}) para el que se desea evaluar el método.
2. Entrenar un modelo de k-means con el conjunto de datos para los números de clústers que van desde 1 hasta k_max. Para cada uno de los modelos se debe obtener la varianza intraclúster (W_k).
3. Generar B veces conjuntos de datos aleatorios, entrenar los mismos modelos de k-means con ellos y obtener la varianza intraclúster promedio (W_k^*).
4. Calcular para cada valor de k la diferencia del logaritmo de la varianza intraclúster gap(k) = \log(W_k^*) - \log(W_k).
5. Obtener la incertidumbre en la estadística de gap (s_k).
6. El valor de k óptimo será aquel que maximice gap(k) y cumpla la condición gap(k) >= gap(k+1) - s_{k+1}.

Implementación del método gap statistics en Python

Una posible implementación del método de gap statistics en Python es la que se muestra a continuación.

import numpy as np
from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist


def gap_statistics(data, k_max, num_sim=10):
    """
    Estima el número óptimo de clusters en un conjunto de datos utilizando el método de Gap Statistics.

    Parámetros:
    -----------
    data : array, shape (n_samples, n_features)
        El conjunto de datos de entrada.
    k_max : int
        El número máximo de clusters a considerar.
    num_sim : int, default=10
        El número de simulaciones aleatorias a utilizar para la distribución de referencia.

    Retorna:
    --------
    k_opt : int
        El número óptimo de clusters para el cojunto de datos
    gap : array, shape (k_max,)
        La estadística de Gap para cada número de clusters k.
    s_k : array, shape (k_max,)
        La desviación estándar de la estadística de Gap para cada número de clusters k.
    """

    # Calcula los rangos de los datos
    ranges = np.apply_along_axis(lambda x: x.max() - x.min(), axis=0, arr=data)

    # Inicializa las matrices para la función de error W_k
    w_k = np.zeros(k_max)
    w_k_star = np.zeros((k_max, num_sim))

    # Calcula la función de error W_k y W_k_star para cada valor de k
    for k in range(1, k_max + 1):
        # Ejecuta el algoritmo de K-means para k clusters
        kmeans = KMeans(n_clusters=k).fit(data)
        # Calcula la función de error W_k para los datos de entrada
        w_k[k - 1] = np.sum(np.min(cdist(data, kmeans.cluster_centers_, 'euclidean'), axis=1)) / data.shape[0]
        # Genera num_sim muestras aleatorias de los datos y calcula la función de error W_k_star para cada muestra
        for sim in range(num_sim):
            random_data = np.random.random_sample(size=data.shape) * ranges + np.min(data, axis=0)
            kmeans = KMeans(n_clusters=k).fit(random_data)
            w_k_star[k - 1, sim] = np.sum(np.min(cdist(random_data, kmeans.cluster_centers_, 'euclidean'), axis=1)) / random_data.shape[0]

    # Calcula la estadística de Gap para cada valor de k
    gap = np.log(w_k_star.mean(axis=1)) - np.log(w_k)

    # Calcula la desviación estándar de la estadística de Gap
    s_k = np.sqrt(1 + 1 / num_sim)*np.std(np.log(w_k_star), axis=1)

    # Encuentra el número óptimo de clusters k
    k_opt = np.argmax(gap) + 1

    # Regresa la estadística de Gap y su desviación estándar
    return k_opt, gap, s_k

En este ejemplo se define la función gap_statistics que implementa el método gap statistics usando las clases de Scikit-learn. La función devuelve el número óptimo de clústeres, los valores Gap y su incertidumbre para cada valor de k. Valores que se pueden usar para representar de forma gráfica los resultados del método.

Ejemplo de uso

El código implementado en la sección anterior se puede evaluar con un conjunto de datos aleatorios para comprobar que el resultado es el esperado.

from sklearn.datasets import make_blobs

# Genera un conjunto de datos de ejemplo con 3 clusters
X, y = make_blobs(n_samples=500, centers=4, n_features=2, random_state=42)

# Calcula la estadística de Gap y su desviación estándar para cada valor de k de 1 a 10
k, _, _ = gap_statistics(X, k_max=10, num_sim=10)

print(k)

4

Pudiendo usar la información que devuelve la función para crear una gráfica.

import matplotlib.pyplot as plt

k_max = 10

k, gap, s_k = gap_statistics(X, k_max=k_max, num_sim=10)

# Grafica la estadística de Gap y su desviación estándar
plt.plot(range(1, k_max + 1), gap, 'bo-', label='Gap')
plt.fill_between(range(1, k_max + 1), gap - s_k, gap + s_k, alpha=0.5)
plt.xlabel('Número de clusters k')
plt.ylabel('Estadística de Gap')
plt.title('Método de Gap Statistics')
plt.legend()
plt.show()

Conclusiones

Uno de los problemas más habituales a la hora de trabajar con k-means, uno de los algoritmos más populares en análisis de clúster, es conocer el número óptimo de clústeres para entrenar el modelo, se puede solucionar mediante el método de gap statistics. Un método que se puede combinar con otro para tener una mejor estimación de cual es los clústeres en los que se debe dividir un conjunto de datos.

Imagen de NakNakNak en Pixabay