A lo largo de esta práctica veremos como aplicar distintas técnicas no supervisadas así como algunas de sus aplicaciones reales:

Clustering clásico: k-means y la regla del codo.

Clustering con formas y feature engineering.
Optimización con reducción de dimensionalidad: t-SNE.
Aplicación: agrupación de documentos.

Para ello vamos a necesitar las siguientes librerías:

import random

import numpy as np
import pandas as pd
from sklearn import cluster      # Algoritmos de clustering.
from sklearn import datasets     # Crear datasets.
from sklearn import manifold     # Algoritmos de reduccion de dimensionalidad.

# Visualizacion.
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline

1. Clustering clásico: k-means y la regla del codo ¶

Vamos a partir de un dataset de clientes en un negocio retail cualquiera (en el fichero pec2_1.p un DataFrame de pandas en formato pickle o pec2_1.csv en formato CSV).

Para cada cliente se cuenta con 3 variables:

n_days_per_week: frecuencia de asistencia a la tienda a la semana.
n_month_purchases: número de compras al mes.
avg_month_turnover: el gasto medio de un cliente al mes.
Primero se pide visualizar las variables para entender como están distribuidas y preprocesarlas para aplicar un k-means.

Implementación: visualizar y preprocesar las variables (para que todas las variables tengan el mismo peso).

df = pd.read_csv("pec2_1.csv")

print("Visualizo las variables y su frecuéncia.  ")
df.hist()
plt.tight_layout()

Visualizo las variables y su frecuéncia.


print("Primer análisis de los datos")
display(df.head())

print("Descripción de los datos")
display(df.describe())
print("El número de líneas es: " + str(df.shape[0]) + " y el número de columnas: "+ str(df.shape[1]))

print("No existe ningún null")
display(df.isnull().sum())

print("Normalizo las variables")
from sklearn.preprocessing import MinMaxScaler
#importante la copia, sino realizo las modificaciones sobre df. 
datos_normalizado = df.copy()
for col in df.columns:
    datos_normalizado[col] = MinMaxScaler().fit_transform( datos_normalizado[col].values.reshape(-1, 1))

display(datos_normalizado.head())

Primer análisis de los datos

Descripción de los datos

El número de líneas es: 1200 y el número de columnas: 3
No existe ningún null

avg_month_turnover    0
n_days_per_week       0
n_month_purchases     0
dtype: int64

Normalizo las variables

Se pide estimar el número de clusters a detectar por k-means. Una técnica para estimar k es, como se explica en la teoría:

Los criterios anteriores (minimización de distancias intra grupo o maximización de distancias inter grupo) pueden usarse para establecer un valor adecuado para el parámetro k. Valores k para los que ya no se consiguen mejoras significativas en la homogeneidad interna de los segmentos o la heterogeneidad entre segmentos distintos, deberían descartarse.

Lo que popularmente se conocer como regla del codo.

Primero es necesario calcular la suma de los errores cuadráticos (SSE) que consiste en la suma de todos los errores (distancia de cada punto a su centroide asignado) al cuadrado.

SSE = \sum_{i=1}^{K} \sum_{x \in C_i} euclidean(x, c_i)^2

Donde K es el número de clusters a buscar por k-means, $x \in C_i$ son los puntos que pertenecen a i-ésimo cluster, $c_i$ es el centroide del cluster $C_i$ (al pertenece el punto x), y euclidean es la distancia euclídea.

Este procedimiento realizado para cada posible valor k, resulta en una función monótona decreciente, donde el eje x representa los distintos valores de k, y el eje y el SSE. Intuitivamente se podrá observar un significativo descenso del error, que indicará el valor idóneo de k.

Se pide realizar la representación gráfica de la regla del codo junto a su interpretación, utilizando la librería matplotlib y la implementación en scikit-learn de k-means.

Implementación: cálculo y visualización de la regla del codo.

# Metodo del Codo para encontrar el numero optimo de clusters

features=['avg_month_turnover', 'n_days_per_week','n_month_purchases']
x = datos_normalizado.loc[:, features].values

from sklearn.cluster import KMeans
wcss = []
x_label = []
for i in range(1, 20):
    kmeans = KMeans(n_clusters = i, init = 'k-means++', random_state = 42)
    kmeans.fit(x)
    wcss.append(kmeans.inertia_)
    x_label.append(i)

# Grafica de la suma de las distancias
plt.plot(range(1, 20), wcss)

plt.title('The Elbow Method')
plt.xlabel('Number of clusters')
plt.xticks(x_label)

plt.ylabel('WCSS')
plt.show()

Análisis: ¿Qué se interpreta en la gráfica? ¿Cómo podría mejorarse la elección de k?.

En la gráfica se observa una disminución muy importante para el valor k=4, es decir, la subdivisión en 4 clústers es la opción más acertada para resolver el problema. A partir de ese valor la disminución del error es muy poco relevante.

Utilizando otra métrica, por ejemplo, se podría utilizar una técnica para maximizar la suma de distancias entre segmentos. Utilizando una tecnica que uniera la minimización de distancias intra- grupo y la maximización de distancias intergrupo obtendríamos mejores resultados.

Análisis: Observando los centroides de cada cluster. ¿Qué tipos de usuarios describen cada cluster?

kmeans = KMeans(n_clusters = 4, init = 'k-means++', random_state = 42)
kmeans.fit(x)
centroids = kmeans.cluster_centers_
print(datos_normalizado.columns)
display(centroids)

Index(['avg_month_turnover', 'n_days_per_week', 'n_month_purchases'], dtype='object')

array([[0.40703499, 0.6464966 , 0.70940991],
       [0.11329995, 0.08032122, 0.05008262],
       [0.64888252, 0.15175766, 0.16899211],
       [0.17023367, 0.51330559, 0.29418709]])

Representan la siguiente tipología de clientes:

El cliente habitual, es decir, realizas compras frecuentes y con un valor elevado.
El cliente ocasionarl, acude muy poco en el mes para comprar y gastar muy poco.
El cliente que gasta mucho, pero que acude muy poco a la tienda.
El cliente que gasta poco, pero acude mucho a la tienda.

[OPCIONAL] Implementación: visualiza el dataset en 3 dimensiones, donde los puntos del mismo color pertenezcan al mismo cluster.

# for 3D projection to work
from mpl_toolkits.mplot3d import Axes3D

import warnings
random.seed(10)
warnings.simplefilter('ignore')

estimators = [('k_means_4', KMeans(n_clusters=4)),
              ]
print('Visualización del dataset en 3 dimensiones')
titles = ['4 clusters']
for name, est in estimators:
    fig = plt.figure().gca(projection='3d')
    est.fit(x)
    labels = est.labels_
    fig.scatter(datos_normalizado['n_days_per_week'], datos_normalizado['n_month_purchases'], datos_normalizado['avg_month_turnover'],
               c=labels.astype(np.float), edgecolor='k')
    fig.set_xlabel('n_days_per_week')
    fig.set_ylabel('n_month_purchases')
    fig.set_zlabel('avg_month_turnover')
    fig.set_title(titles[0])
    plt.tight_layout()
    plt.show()

Visualización del dataset en 3 dimensiones

De forma optativa se plantea realizar el apartado anterior con una implementación propia del algoritmo k-means.

[OPCIONAL] Implementación: algoritmo k-means desde cero.

import copy
import random
random.seed(10)

data=x
# Number of clusters
k = 4
# Number of training data
n = data.shape[0]
# Number of features in the data
c = data.shape[1]

# Generate random centers, here we use sigma and mean to ensure it represent the whole data
mean = np.mean(data, axis = 0)
std = np.std(data, axis = 0)
centers = np.random.randn(k,c)*std + mean



centers_old = np.zeros(centers.shape) # to store old centers
centers_new = copy.deepcopy(centers) # Store new centers

data.shape
clusters = np.zeros(n)
distances = np.zeros((n,k))

error = np.linalg.norm(centers_new - centers_old)

# When, after an update, the estimate of that center stays the same, exit loop
while error != 0:
    # Measure the distance to every center
    for i in range(k):
        distances[:,i] = np.linalg.norm(data - centers[i], axis=1)
    # Assign all training data to closest center
    clusters = np.argmin(distances, axis = 1)

    centers_old = copy.deepcopy(centers_new)
    # Calculate mean for every cluster and update the center
    for i in range(k):
        centers_new[i] = np.mean(data[clusters == i], axis=0)
    error = np.linalg.norm(centers_new - centers_old)
print(centers_new)

[[0.31925093 0.63182734 0.5404736 ]
 [0.20268337 0.27266971 0.2277484 ]
 [0.11899447 0.00887271 0.02935249]
 [0.68292228 0.15964736 0.1608252 ]]

print("Comparación de centroides")
print("Algoritmo kmeans libreria python")
kmeans = KMeans(n_clusters = 4, init = 'k-means++', random_state = 42)
kmeans.fit(x)
centroids = kmeans.cluster_centers_
print(datos_normalizado.columns)
display(centroids)
print("Implementado manualmente")
print(centers_new)

Comparación de centroides
Algoritmo kmeans libreria python
Index(['avg_month_turnover', 'n_days_per_week', 'n_month_purchases'], dtype='object')

array([[0.40703499, 0.6464966 , 0.70940991],
       [0.11329995, 0.08032122, 0.05008262],
       [0.64888252, 0.15175766, 0.16899211],
       [0.17023367, 0.51330559, 0.29418709]])

Implementado manualmente
[[0.31925093 0.63182734 0.5404736 ]
 [0.20268337 0.27266971 0.2277484 ]
 [0.11899447 0.00887271 0.02935249]
 [0.68292228 0.15964736 0.1608252 ]]

print('Visualización de lo centroides')
from mpl_toolkits.mplot3d import Axes3D

titles = ['4 clusters']
for name, est in estimators:
    fig = plt.figure().gca(projection='3d')
    est.fit(x)

    fig.scatter(datos_normalizado['n_days_per_week'], datos_normalizado['n_month_purchases'], datos_normalizado['avg_month_turnover'],
             alpha=.1  )
    fig.scatter(centers_new[:,1], centers_new[:,2],centers_new[:,0], marker='*',alpha=1, c='g', s=[400,400,400,400])

    fig.set_xlabel('n_days_per_week')
    fig.set_ylabel('n_month_purchases')
    fig.set_zlabel('avg_month_turnover')
    fig.set_title(titles[0])
    #for i in range(k):
     #   points = np.array([datos_normalizado[j] for j in range(len(datos_normalizado[j])) if clusters[j] == i])
      #  ax.scatter(points[:, 0], points[:, 1], s=7, c=colors[i])

    plt.tight_layout()
    plt.show()

Visualización de lo centroides

En este caso, con 2 dimensiones, es muy sencillo inferir el número de clusters visualizando los datos. Pero este método es de gran utilidad cuando se cuenta con datos de alta dimensionalidad.

2. Clustering con formas y feature engineering ¶

Pero no todos los datasets son como los del ejercicio anterior. Para esta segunda parte vamos a emplear el siguiente conjunto de datos:

data_circles = ('circles', *datasets.make_circles(n_samples=1000, factor=.5, noise=.05))

Donde data_circles es una tupla con tres posiciones: el nombre del dataset y los dos valores devueltos por la función que genera el dataset:

#datasets.make_circles?

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.scatter(data_circles[1][:,0], data_circles[1][:,1], c=data_circles[2], s=2)
ax.set_title('Dataset {}'.format(data_circles[0]))
plt.tight_layout()

2 a. Encontrando los clusters con k-means¶

Implementación: aplica la regla del codo para decidir el valor de k.

# Metodo del Codo para encontrar el numero optimo de clusters

x=data_circles[1][:, [0,1]]

from sklearn.cluster import KMeans
wcss = []
for i in range(1, 20):
    kmeans = KMeans(n_clusters = i, init = 'k-means++', random_state = 42)
    kmeans.fit(x)
    wcss.append(kmeans.inertia_)
    x_label.append(i)

# Grafica de la suma de las distancias
plt.plot(range(1, 20), wcss)

plt.title('The Elbow Method')
plt.xlabel('Number of clusters')
plt.xticks(range(1, 20))

plt.ylabel('WCSS')
plt.show()

# Metodo del Codo para encontrar el numero optimo de clusters

x=data_circles[1][:, [0,1]]

from sklearn.cluster import KMeans
wcss = []
for i in range(1, 7):
    kmeans = KMeans(n_clusters = i, init = 'k-means++', random_state = 42)
    kmeans.fit(x)
    wcss.append(kmeans.inertia_)
    x_label.append(i)

# Grafica de la suma de las distancias
plt.plot(range(1, 7), wcss)

plt.title('The Elbow Method')
plt.xlabel('Number of clusters')
plt.xticks(range(1, 7))

plt.ylabel('WCSS')
plt.show()

print('La k escogida óptima es 3')

La k escogida óptima es 3

Implementación: aplica k-means con el valor de k elegido.
Visualiza el resultado en un scatter plot representando cada cluster con un color distinto.

# for 3D projection to work
from mpl_toolkits.mplot3d import Axes3D
import warnings
random.seed(10)
warnings.simplefilter('ignore')

estimators = [('k_means_3', KMeans(n_clusters=3)),
              ]
print('Visualización del dataset en 3 dimensiones')
titles = ['3 clusters']
for name, est in estimators:
    fig = plt.figure().gca(projection='3d')
    est.fit(x)
    labels = est.labels_
    fig.scatter( x[:, 1], x[:, 0],
               c=labels.astype(np.float), edgecolor='k')
    fig.set_xlabel('x')
    fig.set_ylabel('y')
    fig.set_title(titles[0])
    plt.tight_layout()
    plt.show()

Visualización del dataset en 3 dimensiones

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

En este punto del ejercicio hemos utilizado Kmeans con distancia euclídea. Esta métrica se basa en la distancia entre los puntos, por lo tanto, el cálculo no tiene en cuenta que existen dos circulos y los tendria que interpretar como diferentes unidades.

Para llevar a cabo este ejercicio sería conveniente utilizar otra métrica, como por densidad para ver si mejora el resultado.

2 b. Más allá de K-Means: algoritmos basados en densidad¶

En este apartado se pide aplicar clustering por densidad como DBSCAN al dataset anterior para poder encontrar los dos clusters iniciales.

Implementación: prueba la implementación de DBSCAN en scikit-learn jugando con los parámetros eps y min_samples para encontrar las 2 estructuras subyacentes (y outliers).

warnings.simplefilter('ignore')

model = cluster.DBSCAN(eps=.1, n_jobs=-1)
fit = model.fit(data_circles[1])

x=data_circles[1][:, [0,1]]
fig = plt.figure(figsize=((10,8))).gca(projection='3d')

labels = fit.labels_
fig.scatter( x[:, 1], x[:, 0],
               c=labels.astype(np.float), edgecolor='k')
fig.set_xlabel('x')
fig.set_ylabel('y')
fig.set_title('DBSCAN')
plt.tight_layout()
plt.show()

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

Con el algoritmo DBSCAN se pueden identificar ambas circunferencias correctamente. Esto es debido a que utilizamos un algoritmo de agrupamiento basado en densidad (density-based clustering) y con los datos proporcionados ha conseguido encontrar un número de grupos (clusters) utilizando la densidad de los puntos.

2 c. Más allá de K-Means: algoritmos jerárquicos¶

En este apartado se pide visualizar mediante un dendrograma la construcción progresiva de los grupos mediante un algoritmo jerárquico aglomerativo (estrategia bottom-up). Con ello se pretende encontrar un método gráfico para entender el comportamiento del algoritmo y encontrar los dos clusters.

Implementación:
prueba la implementación de clustering jerárquico de scipy probando distintos criterios de enlace o linkage permimtiendo identificar los clusters subyacentes (mostrando su resultado) y su dendrograma.

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

fig = plt.figure()
dendrogram(linkage(pdist(x), 'complete') )
plt.title('Denograma con criterios de enlace completo')
plt.ylabel('distance')
plt.tight_layout()
plt.show()

print('La distribución muestras 6 clusters, es decir, no ha funcionado. ')

La distribución muestras 6 clusters, es decir, no ha funcionado.

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

fig = plt.figure()
dendrogram(linkage(pdist(x), 'single') )
plt.title('Denograma con criterios de enlace singl')
plt.ylabel('distance')
plt.tight_layout()
plt.show()

print('La distribución muestras 2 clusters, es decir, lo deseado. ')

La distribución muestras 2 clusters, es decir, lo deseado.

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

fig = plt.figure()
cluster = fcluster(linkage(pdist(data_circles[1]), 'complete'), t=.2, criterion='distance')
plt.scatter( x[:, 1], x[:, 0],
               c=cluster, edgecolor='k')
plt.tight_layout()
plt.show()
print('La distribución muestras 4 clusters. ')

La distribución muestras 4 clusters.

from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster

fig = plt.figure()
cluster = fcluster(linkage(pdist(data_circles[1]), 'single'), t=.2, criterion='distance')
plt.scatter( x[:, 1], x[:, 0],
               c=cluster, edgecolor='k')
plt.tight_layout()
plt.show()
print('Correcto')

Correcto

Análisis: Interpreta el dendrograma y comenta qué criterio de enlace se ha comportado mejor. ¿Por qué?

En nuestro dataset el enlace simple nos ha permitido una mejor clasificación por clusters que el enlace completo. Esto es debido a que el algoritmo ha juntado los puntos más próximos entre sí obteniendo la división deseada. Si los círculos estuvieran más juntos no hubiera sido posible.

Por otra parte, el enlace completo busca máximos, es decir, los que están más alejados, obteniendo como resultado una clasificación peor.

Análisis: ¿Qué ha sucedido? Explica los motivos por los que crees que se ha producido ese resultado.

Ha conseguido encontrar ambos círculos debido a que la densidad es constante y están bastante separados entre ellos.

2 d. Feature engineering y agrupamiento¶

Algunos de los algoritmos anteriores se basan en unas suposiciones que no cumplían en el dataset. Muchas veces en lugar de optar por algoritmos más complejos o que requieren mayor cómputo, se pueden transformar los datos para poder aplicar con éxito técnicas más simples. Esto es un claro ejemplo de feature engineering.

Implementación: transforma los puntos anteriores del dataset a un nuevo espacio de 2 dimensiones:

Radio, o distancia al punto (0,0)
Ángulo,con respecto al vector (1,0)

Para que todas las dimensiones tengan el mismo peso, además vamos a normalizarlas entre 0 y 1 de acuerdo a su máximo y mínimo.
Visualizar los puntos del "nuevo" dataset.

import math

from sklearn.preprocessing import MinMaxScaler

# Distance function
def distance(xi,xii,yi,yii):
       sq1 = (xi-xii)*(xi-xii)
       sq2 = (yi-yii)*(yi-yii)
       return np.sqrt(sq1 + sq2)


radio = np.apply_along_axis(lambda x: distance(x[0], 0,x[1], 0), 1, x)


angulo = np.apply_along_axis(lambda x: math.acos(x[0] / distance(x[0], 0,x[1], 0)), 1, x)

radio = MinMaxScaler().fit_transform(radio.reshape(-1, 1))
angulo = MinMaxScaler().fit_transform(angulo.reshape(-1, 1))

fig = plt.figure()
plt.scatter(radio, angulo, s=4)
plt.tight_layout()
plt.show()

Análisis: ¿Qué crees que sucederá al aplicar los anteriores algoritmos en este "nuevo" dataset?

Ahora si que se puede aplicar la distancia euclídea correctamente, es decir, podemos aplciar el algoritmo sin problemas.

Implementación: aplica cada uno de los algoritmos de agrupamimento anteriores que no hayan podido localizar adecuadamente los dos clusters originales para tratar de encontrarlos en este "nuevo" espacio. Ajusta los parámetros necesarios para facilitar su detección.

Para cada algoritmo, visualiza los clusters encontrados en 2 imágenes:

En el "nuevo" espacio (radio y ángulo).
En el espacio original (posición x e y), pero NO con las etiquetas (pertenencia al cluster) obtenidas al aplicar los algoritmos sobre el dataset original, sino con las etiquetas obtenidas al realizar el clustering en el "nuevo" espacio. A ver si así se consiguen solventar los problemas iniciales.

fig = plt.figure()
from sklearn import cluster      # Algoritmos de clustering.

kmeans = cluster.KMeans(n_clusters=2)
kmeans.fit_transform(np.hstack([radio , angulo]))
labels = kmeans.labels_
plt.scatter(x[:, 1], x[:, 0],
               c=labels.astype(np.float), edgecolor='k')

plt.tight_layout()
plt.show()

3. Aplicación para comprimir imágenes ¶

Las imágenes en color se componen de píxeles que tienen tres componentes (roja, azul y verde), destinando 1 byte a cada canal. Pudiendo representar $2*8 = 256$ rojos, verdes y azules. Con un total de $2^{8^{3}} = 2^{24} = 16777216$ colores representables en cada píxel.

Entre mayor sea el número de colores representables, más memoria será necesaria para almacenar la imagen. Por tanto, una estrategia para comprimir una imagen puede ser disminuir los colores representables en cada píxel, necesitando menos bits para guardar el valor de un píxel. Este método no es reversible, es decir, no se puede volver a recuperar la imagen original a partir de la comprimida. Por tanto, este tipo de compresiónse denomina comprensión con pérdidas.

Pero ¿cómo seleccionamos los "píxeles parecidos" en la imagen original y determinamos su color en la imagen comprimida?. Una opción es utilizar k-means donde k será el número de colores representables, los puntos que pertenecen a cada cluster equivaldrían a los "píxeles parecidos" y las coordenadas de los centroides actuarán como los colores finales a los que se aproximarán los "píxeles parecidos".

Como resultado del clustering, se obtiene una paleta de colores reducida (coordenadas de los centroides) donde cada píxel de la imagen hará referencia a uno de estos colores (cluster al que pertenece). El uso de paletas de colores o colores indexados es un recurso empleado por distintos formatos de imagen como PNG, GIF o TIFF.

Si no dispones de la librería skimage puedes instalarla:

Usando conda (si has creado tu entorno virtual con conda): conda install scikit-image
Usando pip: pip install scikit-image

from skimage import io, transform
import ssl

ssl._create_default_https_context = ssl._create_unverified_context


photo = (transform.resize(
    io.imread('https://upload.wikimedia.org/wikipedia/en/7/7d/Lenna_%28test_image%29.png'),
    (256, 256), mode='edge') * 255).astype(np.uint8)
plt.imshow(photo)        # np.array con shape (256, 256, 3), alto por ancho por 3 (los tres canales: rojo, verde y azul), donde cada valor ocupa un byte.
X = photo.reshape(-1, 3) # np.array con shape (65536, 3), cada pixel con sus 3 canales
print('Imagen con {
                  } pixeles ({} bytes)'.format(X.shape[0], X.shape[0] * 3))

Imagen con 65536 pixeles (196608 bytes)

Para facilitar la tarea, la imagen anterior está guardada en un array con tantas filas como píxeles y columnas como canales (rojo, verde y azul). De forma que cada "muestra" equivale al valor de un píxel.

Se puede volver a recomponer la imagen original con:

plt.imshow(X.reshape(photo.shape))

Podemos ver cada píxel como un punto en un sistema de coordenadas de 3 dimensiones donde una es su cantidad de rojo, otra su cantidad de verde y por último su cantidad de azul. Por lo que podemos realizar una visualización en 3 dimensiones de los píxeles sobre la que se probará el algoritmo de clustering:

X

array([[224, 136, 125],
       [223, 136, 128],
       [224, 132, 118],
       ...,
       [165,  68,  82],
       [175,  67,  78],
       [181,  70,  79]], dtype=uint8)

Para visualizar la animación de la representación 3D de los píxeles es necesario instalar la librería ffmpeg.

Si tienes un entorno virtual de conda, lo puedes instalar con: conda install -c conda-forge ffmpeg.

Si tienes un error, puedes ver una representación estática 2 celdas más adelante.

Visualización estática del valor de los píxeles:

idx = np.random.randint(X.shape[0], size=int(X.shape[0] * .3))
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig)
ax.scatter(X[idx,0], X[idx,1], X[idx,2], s=1, c=X[idx,:] / 255, alpha=.3)
ax.set_xlabel('rojo')
ax.set_ylabel('verde')
ax.set_zlabel('azul')
ax.view_init(20, 350)
plt.show()

En la visualización anterior se ha representado cada píxel con su color, donde sus coordenadas en los 3 colores oscilan entre 0 (carece de esa componente) y 1. Podemos comprobar como los píxeles en coordenadas (255, 255, 255) son píxeles blancos y los situados en (0, 0, 0) son píxeles negros. También se observan 4 estructuras de píxeles: la pluma del sombrero, el sombrero, la piel de la modelo (Lenna) y el fondo.

Implementación: Aplicar k-means con k=3 a los píxeles de la imagen (puntos con 3 dimensiones de la variable X) y obtener, para cada punto, su centroide más cercano. Y las coordenadas (3 dimensiones) de cada centroide.

from PIL import Image
from sklearn.cluster import KMeans
im = Image.open("lena.png")
# display image
display(im)
# get pixels of the image
pixel_np = np.asarray(im)
# reshape array (remove rows and columns)
image_height = im.height
image_width = im.width
pixel_np = np.reshape(pixel_np, (image_height * image_width, 3))
# display as df
display(pd.DataFrame(pixel_np, columns=["r", "g", "b"]).head())

# run k-means clustering on the pixel data
num_of_centroids = 3#3 # a 4-bit image is represented by 2^4 colours
num_of_runs = 10 # number of times to run the k-means algorithm before determining the best centroids
max_iterations = 300 # number of iterations before k-means comes to an end for a single run
verbosity = 0 # show what's going on when the algorithm is running

# initiate a kmeans object
compressor = KMeans(n_clusters=num_of_centroids, n_init=num_of_runs, max_iter=max_iterations, verbose=verbosity)
# run k-means clustering
compressor.fit(pixel_np)

display(compressor.cluster_centers_)

array([[224.83959878, 166.770715  , 147.58069924],
       [110.05040217,  34.4708127 ,  71.89841545],
       [195.41033953,  97.96534435, 100.87512045]])

Implementación: asgina a cada punto (píxel de la imagen) el valor de su centroide asociado en lugar de su propio valor.

# create an array replacing each pixel label with its corresponding cluster centroid
pixel_centroid = np.array([list(compressor.cluster_centers_[label]) for label in compressor.labels_])

Implementación: volver a convertir el resultado en una imagen usando la función reshape (como cuando se creó X) para que vuelva a tener su dimensión de 256 x 256. Y mostrarla con imshow.

# convert the array to an unsigned integer type
pixel_centroid = pixel_centroid.astype("uint8")
# reshape this array according to the height and width of our image
pixel_centroids_reshaped = np.reshape(pixel_centroid, (image_height, image_width, 3), "C")




# create the compressed image
compressed_im = Image.fromarray(pixel_centroids_reshaped)
# save compressed image
compressed_im.save("bridge_compressed.jpeg")

plt.figure()
plt.imshow(pixel_centroids_reshaped)
plt.axis('off')
plt.title('3 clusters')
plt.tight_layout()
plt.show()

Análisis: ¿por qué se ha producido este resultado? ¿qué relación tiene con el clústering?

Por que hemos aplicado una cuantificación, es decir, hemos aplicado una técnica de compresión con pérdida que consiste en agrupar todo un rango de valores en uno solo. Si cuantificamos el color de una imagen, reducimos el número de colores necesarios para representarla y el tamaño del fichero de la misma disminuye.

El clústering es la técnica que nos permite llevar a cabo la compresión.

Implementación: realizar el proceso anterior para distintos valores de k: 256, 128, 64, 32, 16, 8, 4 y 2 colores. Mostrar las imágenes e indicar los bytes que ocuparía cada una de las opciones.

from PIL import Image


for clusters_utilizar in  [ 128, 64, 32, 16, 8, 4 , 2,256]:
    display(clusters_utilizar)
    im = Image.open("lena.png")
    # get pixels of the image
    pixel_np = np.asarray(im)
    # reshape array (remove rows and columns)
    image_height = im.height
    image_width = im.width
    pixel_np = np.reshape(pixel_np, (image_height * image_width, 3))
    # display as df
    # run k-means clustering on the pixel data
    num_of_centroids = clusters_utilizar
    num_of_runs = 10 # number of times to run the k-means algorithm before determining the best centroids
    max_iterations = 300 # number of iterations before k-means comes to an end for a single run
    verbosity = 0 # show what's going on when the algorithm is running

    # initiate a kmeans object
    compressor = KMeans(n_clusters=num_of_centroids, n_init=num_of_runs, max_iter=max_iterations, verbose=verbosity)
    # run k-means clustering
    compressor.fit(pixel_np)

    display(compressor.cluster_centers_)
    # convert the array to an unsigned integer type
    pixel_centroid = pixel_centroid.astype("uint8")
    # reshape this array according to the height and width of our image
    pixel_centroids_reshaped = np.reshape(pixel_centroid, (image_height, image_width, 3), "C")




    # create the compressed image
    compressed_im = Image.fromarray(pixel_centroids_reshaped)

    plt.figure()
    plt.imshow(pixel_centroids_reshaped)
    plt.axis('off')
    plt.title(clusters_utilizar)
    plt.tight_layout()
    display(plt.show())

128

array([[ 87.73353846,  20.21353846,  63.46482051],
       [210.37828869, 137.0123252 , 122.8113297 ],
       [180.45664967,  78.14123186,  89.01137701],
       [228.47013889, 205.39444444, 194.59375   ],
       [216.82851344, 111.77067478,  96.47590056],
       [100.41058394,  36.77919708,  87.82360097],
       [235.01756757, 149.46722973, 121.17635135],
       [137.0437018 ,  86.697515  , 139.77806341],
       [121.44133644,  52.29292929,  92.92851593],
       [221.91274298, 183.30583153, 166.75637149],
       [219.33767442, 105.9655814 , 111.61534884],
       [103.83255471,  29.25178264,  65.28964839],
       [161.61680551,  54.47229265,  71.81492213],
       [159.38807786, 127.25912409, 170.79075426],
       [239.60292524, 201.20693391, 157.74702059],
       [202.93303065,  95.60017026,  98.77525539],
       [226.10088272, 126.03656999, 103.04319042],
       [209.35343788, 146.09770808, 146.38962606],
       [185.48101983, 111.11558074, 111.25269122],
       [124.42001576,  64.84239559, 117.85657998],
       [194.37114498, 127.87344913, 124.71960298],
       [140.84596211,  40.26570289,  65.40877368],
       [152.90007686,  77.18831668, 100.45810915],
       [235.62654668, 148.39932508, 109.68841395],
       [174.65972222, 110.71219136, 119.29089506],
       [190.813273  ,  80.55656109,  86.08295626],
       [240.93874834, 182.94141145, 121.7310253 ],
       [173.9202454 ,  66.39306114,  78.46498836],
       [213.93642827,  89.48347972,  85.98494354],
       [193.81398417, 160.82717678, 162.3469657 ],
       [231.15719523, 135.69019248, 103.80018332],
       [ 94.13501787,  27.53726391,  69.89739663],
       [114.33981002,  36.4343696 ,  68.50086356],
       [232.53447185, 177.15686275, 157.14611006],
       [225.40115942, 116.10318841,  94.80289855],
       [237.6051109 , 158.81533269, 127.59980714],
       [241.26537217, 206.77605178, 188.98446602],
       [221.672748  , 105.96579247, 101.30359179],
       [243.69173729, 191.82733051, 135.81673729],
       [174.92239468, 102.71175166, 139.63192905],
       [224.55855397, 127.27393075, 125.72454175],
       [194.98291925, 121.9673913 , 114.64518634],
       [225.4834653 , 192.07591989, 176.870517  ],
       [171.75126904, 101.70642978, 106.51692047],
       [212.73819055, 158.66933547, 159.58366693],
       [185.22733954,  86.2199902 ,  95.93483586],
       [107.73106865,  45.69072895,  97.15145081],
       [ 92.62301396,  30.73952817,  79.57775638],
       [146.70548367,  67.978435  ,  88.12322859],
       [208.61267981, 112.91209377, 114.10069259],
       [248.08115942, 208.84492754, 156.61666667],
       [152.30726016,  47.4187554 ,  68.80812446],
       [213.43469786, 173.28330084, 169.02469136],
       [202.59246213, 134.83127862, 127.47516731],
       [221.5778412 , 116.2540218 , 106.51842242],
       [116.46545455,  26.96606061,  59.59878788],
       [238.86977887, 218.54668305, 203.15724816],
       [159.98913043,  73.25      ,  87.6297954 ],
       [193.27837965,  74.17475271,  77.75459256],
       [231.37162162, 154.04560811, 145.8277027 ],
       [ 93.29479769,  16.79710983,  56.97023121],
       [239.11888587, 205.92527174, 169.58152174],
       [216.24484305,  94.46726457,  99.90941704],
       [ 73.04985673,   9.54383954,  55.02063037],
       [191.26907631,  95.01204819, 100.39848282],
       [195.91457286, 127.76884422, 143.52135678],
       [172.90963257,  71.68520357,  85.75571003],
       [226.21619667, 168.80187997, 144.4692697 ],
       [182.15559157, 143.39546191, 151.35656402],
       [154.14383562,  98.18835616, 129.81506849],
       [127.17664975,  45.69847716,  72.70050761],
       [222.24675325, 155.65481887, 134.04511278],
       [151.83786317, 112.84629803, 160.96157451],
       [134.79615385,  64.36730769, 100.88942308],
       [224.88572933, 138.03159558, 122.63717746],
       [236.29651861, 186.47839136, 169.38895558],
       [195.93642305,  88.69073084,  91.16666667],
       [168.19906323,  83.77751756,  96.75800156],
       [102.65048895,  34.65121333,  75.14270192],
       [180.85598142,  90.45412311, 114.76422764],
       [202.41832267,  83.10811721,  82.78275514],
       [136.19871383,  57.72990354,  81.42893891],
       [ 78.2578304 ,  16.22383499,  64.73834989],
       [209.03995902, 122.37397541, 102.01844262],
       [164.03456298,  62.97308076,  79.29145896],
       [210.3356586 , 137.01148954, 134.00533443],
       [219.9331761 , 101.99095912,  90.81092767],
       [217.93079179, 185.31612903, 178.42991202],
       [234.86120401, 181.04598662, 133.45819398],
       [208.40756994, 123.2024136 , 125.15085025],
       [214.42611684, 159.91271478, 147.41030928],
       [144.54390452,  99.45609548, 151.27109974],
       [103.79027903,  21.1039604 ,  57.54365437],
       [128.35910364,  35.19831933,  63.31204482],
       [195.09407666, 110.90679443, 123.88501742],
       [223.85463594, 129.34101624, 112.39130435],
       [205.46003717,  90.99790892,  90.17611524],
       [236.42183463, 160.07364341, 117.40762274],
       [ 84.23982737,  13.39303329,  56.0782984 ],
       [160.28976234,  89.64351005, 110.91316271],
       [231.85984696, 139.46999597, 113.92871526],
       [ 96.40490798,  23.31792287,  62.80324277],
       [140.38039944,  49.80817464,  72.51881096],
       [170.75478261, 123.33217391, 138.30086957],
       [244.46721311, 199.40710383, 146.83196721],
       [185.34737874, 120.77756002, 121.44439   ],
       [151.85976409,  58.600699  ,  77.42027086],
       [212.31989833, 145.42193174, 129.32207698],
       [115.74862096,  55.15681639, 109.10480693],
       [139.58867925,  79.32075472, 119.14968553],
       [249.55936675, 216.38258575, 169.27968338],
       [ 85.03164557,  23.23194341,  72.11057334],
       [181.10767947,  64.95868114,  73.81051753],
       [167.86761711, 145.80448065, 183.92260692],
       [187.92857143, 172.75824176, 193.54945055],
       [208.15161204,  98.81188383,  92.94191314],
       [238.78827362, 197.20651466, 181.12899023],
       [204.44115031, 130.7269493 , 117.90661133],
       [128.9540567 ,  73.62756598, 131.40762463],
       [227.3952591 , 141.64167585, 135.53693495],
       [224.77671451, 196.70095694, 186.33253589],
       [207.63098986, 103.82896118, 103.22455404],
       [217.44053601, 172.1345617 , 157.07481854],
       [221.56787834, 118.57900593, 117.34681009],
       [171.96600385,  60.02148813,  72.50801796],
       [198.20823389, 101.69212411, 110.52505967],
       [181.87722222,  72.15527778,  80.72694444],
       [111.52333136,  43.37684584,  82.87714117]])

None

64

array([[229.85576317, 135.18769683, 106.85387361],
       [147.76064356,  46.58094059,  69.04826733],
       [193.61351882,  87.0062422 ,  91.05243446],
       [243.85789938, 208.74786702, 164.48426008],
       [ 84.79679767,  20.51033479,  66.52372635],
       [211.07449179, 140.60818713, 128.13659148],
       [141.14142539,  91.2811804 , 141.42928731],
       [220.79424407, 111.62549056,  95.68155485],
       [228.81221198, 176.0437788 , 154.97119816],
       [165.96708464,  68.78526646,  84.00515002],
       [102.13263212,  37.64674192,  84.63442791],
       [207.31767417,  87.29161496,  85.53837972],
       [202.0477126 , 112.99719338, 115.04911591],
       [113.45842843,  26.77698076,  59.92957287],
       [234.875     , 187.64264706, 169.25294118],
       [131.97613365,  75.28210024, 127.19761337],
       [225.5829447 , 161.        , 142.29913391],
       [159.2092415 , 126.63993025, 170.31211857],
       [219.76524834, 108.03753438, 106.87849862],
       [159.17878513,  56.21795104,  74.02484134],
       [238.95800317, 184.42908082, 130.54160063],
       [210.97312224, 147.42120766, 144.6185567 ],
       [224.56955435, 131.49833403, 121.34089963],
       [191.97073519, 123.78283369, 119.96556031],
       [116.88837517,  39.63606341,  71.68064729],
       [178.90816327,  79.22576531,  90.6860119 ],
       [237.21586716, 214.46432964, 199.38376384],
       [188.85897436,  94.18478261, 103.4099777 ],
       [165.27122836,  93.92085738, 112.17436109],
       [113.74722766,  49.03979126,  95.66699282],
       [144.31591963,  64.71344668,  86.95177743],
       [227.20961659, 144.21998743, 134.52796983],
       [149.10025543, 107.37803321, 157.12452107],
       [ 76.28804348,  11.71813241,  57.45998024],
       [179.42474656,  69.78281778,  79.63737978],
       [205.20057307,  95.25351937,  93.86508035],
       [236.6740757 , 158.11245599, 123.19828345],
       [205.5147257 , 101.59230029, 103.46525505],
       [173.03766478, 153.48399247, 187.00564972],
       [223.66333073, 123.77107115, 104.90742768],
       [243.65356711, 200.61003628, 149.93258767],
       [181.3824228 , 110.57957245, 115.23145949],
       [173.72146447, 115.51399856, 136.74012922],
       [205.77673496, 131.75611368, 120.46609385],
       [172.96135748,  62.13874234,  74.42877513],
       [199.89467312, 130.34231235, 132.94279661],
       [ 97.4805382 ,  22.88779433,  61.88971648],
       [225.60655738, 198.47851894, 187.87478802],
       [234.41073961, 146.80952381, 115.22208713],
       [192.97994601,  77.15387582,  81.13420748],
       [212.11143505, 163.4829932 , 157.48202138],
       [217.69370892,  97.61915493,  93.89652582],
       [122.84953233,  61.15209435, 111.57055714],
       [103.12585278,  31.15888689,  68.46983842],
       [154.74169884,  78.49305019,  99.23204633],
       [187.35269001, 148.59521776, 155.58923997],
       [ 89.7133853 ,  15.85801957,  57.21164372],
       [132.5282838 ,  37.79482263,  64.55608821],
       [133.79228395,  51.56203704,  76.39197531],
       [220.14580941, 118.39425947, 117.29919633],
       [221.47202073, 188.61243523, 177.08445596],
       [239.61745429, 202.25388065, 183.81752328],
       [216.65348007, 177.42028554, 167.57555027],
       [ 91.98831291,  28.24655985,  74.35777568]])

None

32

array([[220.27535507, 119.03660732, 114.48539708],
       [ 95.058863  ,  29.83002588,  74.30227207],
       [223.19205027, 177.87091927, 163.14642808],
       [139.85714286,  47.57674471,  71.44734759],
       [193.45095681,  79.29666484,  83.23247676],
       [213.2906525 , 139.06770792, 129.55284231],
       [151.28717201, 107.27368805, 154.08965015],
       [ 97.92128386,  22.80771186,  61.58158449],
       [136.01805647,  82.05121471, 132.9425476 ],
       [185.59897781, 115.47606582, 119.32136624],
       [206.54216305, 158.03550572, 156.37187573],
       [162.02358143,  57.43509843,  74.12169702],
       [235.40538453, 152.47191879, 120.09875317],
       [225.20010532, 191.52395998, 179.71814113],
       [204.67537443, 102.43898758, 105.09404389],
       [167.03769841,  93.12921627, 109.82142857],
       [208.35379826,  92.86946155,  90.89396649],
       [227.88903096, 132.81295095, 108.59548817],
       [188.97562893,  89.04363208,  96.1889413 ],
       [225.10959154, 156.75929978, 141.75820569],
       [123.03057603,  60.7035217 , 110.16926017],
       [117.47175379,  34.10876897,  65.60160202],
       [168.4368071 , 137.33370288, 170.56152993],
       [201.19227406, 129.12781178, 122.55321517],
       [ 82.68348452,  15.86820601,  60.69487968],
       [239.51367245, 185.76653925, 133.19553073],
       [109.50176482,  43.90469998,  89.27995542],
       [220.75859883, 108.84537962,  99.11911555],
       [236.3580176 , 206.93608879, 191.11174895],
       [176.0821599 ,  68.83406921,  80.23383055],
       [150.64790024,  70.69759765,  91.32000734],
       [243.28277946, 204.55649547, 159.43746224]])

None

16

4. Aplicación de reducción de dimensionalidad para resolver un problema de optimización: t-SNE ¶

Como ya se ha visto, el algoritmo t-SNE ideado por van der Maaten y Hinton difiere de PCA en que no trata de maximizar la varianza explicada. Intuitivamente, t-SNE trata de que la vecindad de un punto en baja dimensionalidad sea la misma que la original (mantenga las distancias). Partiendo de una localización aleatoria de cada punto, corrige su posición de forma iterativa tratando de minimizar la distancia a sus vecinos originales hasta converger.

Para ello, t-SNE dispone de diversos parámetros que pueden modificar drásticamente el resultado. Por lo que se recomienda conocer su funcionamiento antes de aplicar la técnica.

Partiendo de las distancias entre las provincias de la península ibérica, presentes en el fichero de datos (en pec2_2.p un DataFrame de pandas en formato pickle o pec2_2.csv en formato CSV). Se pide calcular la matriz cuadrada que contenga la distancia de cada provincia contra las demás.

Implementación: la matriz debe tener tantas filas y columnas como provincias. Y cada celda debe contener la distancia entre las provincias de esa fila y columna.

df = pd.read_pickle('pec2_2.p')

provincias = np.unique(df['from'].values).tolist()
distancia = np.zeros((len(provincias), len(provincias)))

df.set_index(['from', 'to'], inplace=True)
for i, a in enumerate(provincias):
    for j, b in enumerate(provincias):
        distancia[i, j] = df.loc[a].loc[b].values[0]

Una vez que se cuenta con la matriz de distancias, t-SNE tratará de mantener esas distancias entre los distintos puntos en baja dimensionalidad (en este caso 2 dimensiones). Emplazando los puntos en el plano mientras intenta mantener las distancias indicadas.

Dado que la entrada a t-SNE se le pasa la matriz de distancias, no es necesario que las calcule. Por ello le indicaremos que la métrica a emplear es "precalculada".

Como t-SNE es un algoritmo estocástico (dos ejecuciones consecutivas con los mismos datos pueden conducir a resultados diferentes). Se pide realizar el proceso de ajuste con t-SNE 100 veces y quedarse con la ejecución con menor error (ver el atributo kl_divergence).

Implementación: jugar con los hiperparámetros de t-SNE, ejecutar 100 veces el ajuste de t-SNE y guardar el resultado de emplazamiento de provincias en el plano de la ejecución con menor error.

min_error = 1
coords = None
for i in range(100):
    tsne = manifold.TSNE(n_components=2, metric='precomputed', perplexity=30, n_iter=10000)
    new_coord = tsne.fit_transform(np.array(distancia))
    if tsne.kl_divergence_ < min_error:
        coords = new_coord
        min_error = tsne.kl_divergence_

Una vez que se tienen las posiciones de las provincias en el plano, visualizar el resultado y analizar si el emplazamiento de las provincias calculado por t-SNE calculado en base a las distancias se parece al real.

Implementación: visualizar en un scatter las provincias de la ejecución con menor error, junto con su nombre para poder analizarlo.

fig, ax = plt.subplots()
ax.scatter(coords[:,0], coords[:,1])
for i, ciudad in enumerate(provincias):
    ax.annotate(ciudad, (coords[i,0], coords[i,1]))
plt.show()

Análisis: ¿se parece a la distribución de provincias real? ¿por qué?

Sí, el algoritmo ha mantenido las distancias en el mapa. Porque sabe interpretar la distancia entre las provincia, aunque no sabe como están colocadas realmente en el mapa, es decir, aparece rotado con respecto a la realidad.

	avg_month_turnover	n_days_per_week	n_month_purchases
0	108.331434	3.263182	5.523083
1	101.253142	5.197333	5.303098
2	60.900002	3.176871	5.082450
3	90.926811	4.184747	5.272594
4	114.936497	3.335393	4.754987

	avg_month_turnover	n_days_per_week	n_month_purchases
count	1200.000000	1200.000000	1200.000000
mean	215.487550	2.711946	5.120623
std	180.534793	1.969671	4.362437
min	-41.444818	0.100000	0.100000
25%	79.262486	0.854183	1.578684
50%	133.655287	2.728147	4.233838
75%	331.432495	4.331359	6.857876
max	793.734593	7.685392	16.856447

	avg_month_turnover	n_days_per_week	n_month_purchases
0	0.179334	0.417010	0.323642
1	0.170859	0.671993	0.310513
2	0.122542	0.405631	0.297345
3	0.158495	0.538502	0.308693
4	0.187243	0.426529	0.277803

Clustering y aprendizaje no supervisado: k-means, t-SNE y documentos

Respuesta rápida

k-means

Regla del codo

t-SNE

Documentos e imágenes

Conclusión práctica

Métodos no supervisados

1. Clustering clásico: k-means y la regla del codo ¶

2. Clustering con formas y feature engineering ¶

2 a. Encontrando los clusters con k-means¶

2 b. Más allá de K-Means: algoritmos basados en densidad¶

2 c. Más allá de K-Means: algoritmos jerárquicos¶

2 d. Feature engineering y agrupamiento¶

3. Aplicación para comprimir imágenes ¶

4. Aplicación de reducción de dimensionalidad para resolver un problema de optimización: t-SNE ¶

	r	g	b
0	226	137	125
1	226	137	125
2	223	137	133
3	223	136	128
4	226	138	120