Práctica 1: Algoritmo SOM de Kohonen en una dimension

Los mapas auto-organizativos originalmente propuestos por Teuvo Kohonen se implementaron sobre redes bidimensionales. De hecho, son los mapas más utilizados porque esta dimensión permite una visualización adecuada de los clusters con un esfuerzo computacional razonable.

No obstante, los fundamentos del algoritmo de Kohonen se pueden comprender de una manera si cabe más inmediata si las neuronas forman una red unidimensional. En esta primera sección vamos a aplicar un SOM unidimesional para clasificar una colección de vectores binarios.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Definición de los datos

Definimos una colección de \(n_i\) vectores de entrada (instancias) de dimensión \(n_d\) formados por elementos en el conjunto \(\{0,1\}\). Un parámetro que podemos dejar libre es la proporción de cada uno de los dígitos. Usaremos una proporción de 5 a 1.

elementos = [0,1]
proporcion = [5/6,1/6]

np.random.choice(elementos,size=10,replace=True,p=proporcion)

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1])

# generamos datos iniciales

nd = 10
ni = 1000
columnas = [i for i in np.arange(nd)]
datosIniciales = pd.DataFrame(0,columns = columnas,index=np.arange(ni))
for i in np.arange(ni):
    datosIniciales.loc[i] = np.random.choice(elementos,size=10,replace=True,p=proporcion)
datosIniciales

	0	1	2	3	4	5	6	7	8	9
0	0	0	0	0	0	0	0	0	0	1
1	0	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0	0
3	1	0	0	0	0	0	0	1	1	1
4	0	0	0	0	0	1	0	0	0	0
...	...	...	...	...	...	...	...	...	...	...
995	0	0	0	0	1	0	0	0	0	0
996	0	0	0	0	0	0	0	1	0	0
997	0	0	0	0	0	1	0	0	0	0
998	0	0	0	0	0	0	0	0	1	0
999	1	0	0	0	0	0	0	0	0	0

1000 rows × 10 columns

# sumamos por filas
datosIniciales.sum(axis=1)

0      1
1      0
2      0
3      4
4      1
      ..
995    1
996    1
997    1
998    1
999    1
Length: 1000, dtype: int64

datosIniciales.sum(axis=1).mean()

1.668

datosIniciales.sum().sum(),datosIniciales.sum().sum()/(datosIniciales.shape[0]*datosIniciales.shape[1])

(1668, 0.1668)

1/6

0.16666666666666666

datosIniciales.mean(),datosIniciales.std()

(0    0.171
 1    0.167
 2    0.158
 3    0.178
 4    0.168
 5    0.161
 6    0.160
 7    0.167
 8    0.169
 9    0.169
 dtype: float64,
 0    0.376697
 1    0.373162
 2    0.364924
 3    0.382704
 4    0.374053
 5    0.367715
 6    0.366789
 7    0.373162
 8    0.374939
 9    0.374939
 dtype: float64)

f,ax=plt.subplots(nrows=1,ncols=2,figsize=(10,5))
allDataVector = datosIniciales.to_numpy().ravel()
ax[0].hist(allDataVector,bins=20)
ax[0].set_title('Todos los datos iniciales')
ax[1].hist(datosIniciales.loc[0])
ax[1].set_title('Primer dato inicial')
plt.show()

f,ax=plt.subplots(nrows=1,ncols=1,figsize=(5,5))
datosIniciales[:10].plot.bar(stacked=True,ax=ax)

<Axes: >

datosIniciales[:10]

	0	1	2	3	4	5	6	7	8	9
0	0	0	0	0	0	0	0	0	0	1
1	0	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0	0
3	1	0	0	0	0	0	0	1	1	1
4	0	0	0	0	0	1	0	0	0	0
5	0	0	0	0	0	1	0	0	0	0
6	0	0	1	0	0	0	0	0	0	0
7	0	0	0	1	1	0	0	0	0	0
8	0	0	0	0	0	0	0	0	0	0
9	0	0	0	0	1	0	0	0	0	0

f,ax=plt.subplots(nrows=1,ncols=1,figsize=(30,15))
datosIniciales.plot.bar(stacked=True,ax=ax,legend=False,xticks=[])

<Axes: >

datosIniciales.sum(axis=1).max()

6

datosIniciales.sum(axis=1).argmax()

42

Medias y varianzas de los datos

• Por atributo

• Por dato

f,ax=plt.subplots(nrows=1,ncols=2,figsize=(10,5))
x=np.arange(nd)
ax[0].errorbar(x,datosIniciales.mean(),datosIniciales.std(),linestyle='None', marker='^')
ax[0].grid()
ax[0].set_title('Por atributo')
x = np.arange(ni)
ax[1].errorbar(x,datosIniciales.mean(axis=1),datosIniciales.std(axis=1),linestyle='None', marker='^',alpha=0.1)
ax[1].grid()
ax[1].set_title('Por dato')
f.suptitle('Medias y desviaciones típicas')

Text(0.5, 0.98, 'Medias y desviaciones típicas')

datosIniciales.std(axis=1).min(),datosIniciales.std(axis=1).max()

(0.0, 0.5270462766947299)

Inicialización del SOM

En el siguiente paso definimos la red neuronal (número de neuronas, \(neu\)) y asignamos valores a los vectores de referencia asociados a cada neurona, con igual dimensión que las instancias, \(n_d\).

Esta asignación se realiza aleatoriamente y, en la medida de lo posible, debería tener en cuenta la forma de los datos que se van a analizar.

neu = 25

W = np.random.uniform(low=0.0,high=1.0,size=(neu,nd))
Wdf = pd.DataFrame(W)
Wdf.head()

	0	1	2	3	4	5	6	7	8	9
0	0.533451	0.070169	0.422484	0.005212	0.169474	0.500415	0.496419	0.735987	0.604227	0.204534
1	0.311411	0.622224	0.427547	0.878344	0.854660	0.387227	0.630206	0.990831	0.934898	0.421978
2	0.616974	0.695730	0.892648	0.050825	0.555735	0.921758	0.303803	0.816286	0.032678	0.804537
3	0.823512	0.979943	0.331533	0.931243	0.850920	0.215364	0.867593	0.742559	0.608166	0.978343
4	0.176341	0.307356	0.533940	0.493924	0.220288	0.300183	0.385594	0.599003	0.167756	0.161135

Cuarto atributo

Wdf.loc[:,4]

0     0.169474
1     0.854660
2     0.555735
3     0.850920
4     0.220288
5     0.453909
6     0.551690
7     0.356566
8     0.038156
9     0.346988
10    0.951792
11    0.915266
12    0.231304
13    0.707766
14    0.473758
15    0.896741
16    0.334188
17    0.888724
18    0.226369
19    0.062289
20    0.178548
21    0.894781
22    0.460768
23    0.285485
24    0.017471
Name: 4, dtype: float64

Cuarta neurona

Wdf.loc[4,:]

0    0.176341
1    0.307356
2    0.533940
3    0.493924
4    0.220288
5    0.300183
6    0.385594
7    0.599003
8    0.167756
9    0.161135
Name: 4, dtype: float64

Tenemos 10 features/atributos. En los datos éstas toman valores en el conjunto \(\{0,1\}\), mientras que en las neuronas, los pesos toman valores en el intervalo \([0,1]\). Hacemos un histograma de cómo se distribuyen tanto en los datos como en la red neuronal. Tomamos un feature al azar.

Wdf.shape

(25, 10)

Histogramas

• Todos los valores

• Un feature aleatorio

• Un dato vs una neurona aleatorios

f,ax=plt.subplots(nrows=1,ncols=2,figsize=(10,5))
allDataVector = datosIniciales.to_numpy().ravel()
ax[0].hist(allDataVector,bins=20,range=[0,1])
ax[0].set_title('Datos iniciales')
allWVector = Wdf.to_numpy().ravel()
ax[1].hist(allWVector,bins=40,range=[0,1])
ax[1].set_title('SOM inicial')
f.suptitle('Histogramas')
plt.show()

allWVector[:5]

array([0.53345141, 0.07016899, 0.42248415, 0.00521206, 0.1694736 ])

Histograma por feature/atributo

feature = np.random.choice(nd)
datosInicialesFeature = datosIniciales.loc[:,feature]
datosNeuronaFeature = Wdf.loc[:,feature]
f,ax = plt.subplots(nrows=1,ncols=2,figsize=(10,5))
datosInicialesFeature.hist(ax=ax[0],range=[0,1])
datosNeuronaFeature.hist(ax=ax[1],range=[0,1])
ax[0].set_title('Datos iniciales')
ax[1].set_title('SOM inicial')
plt.show()
print('Feature: ',feature)

Feature:  3

Histograma de dato/neurona individual

neurona = np.random.choice(neu)
dato = np.random.choice(ni)
datosInicialesDato = datosIniciales.loc[dato,:]
datosNeurona = Wdf.loc[neurona,:]
f,ax = plt.subplots(nrows=1,ncols=2,figsize=(10,5))
datosInicialesDato.hist(ax=ax[0],range=[0,1])
datosNeurona.hist(ax=ax[1],range=[0,1])
ax[0].set_title('Dato')
ax[1].set_title('Neurona')
plt.show()
print('Dato: {}\nNeurona: {}'.format(dato,neurona))

Dato: 463
Neurona: 9

Medias y varianzas de las neuronas

f,ax=plt.subplots(nrows=1,ncols=2,figsize=(10,5))
x=np.arange(nd)
ax[0].errorbar(x,Wdf.mean(),Wdf.std(),linestyle='None', marker='^')
ax[0].grid()
x = np.arange(neu)
ax[1].errorbar(x,Wdf.mean(axis=1),Wdf.std(axis=1),linestyle='None', marker='^',alpha=0.1)
ax[1].grid()
ax[0].set_title('Por atributo')
ax[1].set_title('Por neurona')

Text(0.5, 1.0, 'Por neurona')

Stacked barplot de las neuronas del SOM en su estado inicial

f,ax=plt.subplots(nrows=1,ncols=1,figsize=(30,15))
Wdf.plot.bar(stacked=True,ax=ax,legend=False)

<Axes: >

Funciones de núcleo de entorno

En la actualización de los vectores de referencia de las neuronas se utiliza un núcleo de entorno (neighborhood kernel), \(h_{ci}\), que modifica el entorno de vectores cercanos a la BMU. En este núcleo aparecen las funciones de la tasa de aprendizaje \(\alpha\) y de varianza \(\sigma\). Una sencilla elección es la siguiente:

\[ \alpha(k)=\frac{1}{1+k}\qquad \sigma(k)=-e^{-0.3k} \]

# graficamos las elecciones
kk = np.arange(0,100)
def alpha(k):
    s = k/(neu)
    return 1/(1+s)
def sigma(k):
    s=k/neu
    return np.exp(-0.3*s)
f,ax=plt.subplots(nrows=1,ncols=2,figsize=(10,5))
ax[0].plot(kk,alpha(kk),'--',label=r'$\alpha(k)$')
ax[1].plot(kk,sigma(kk),'--',label=r'$\sigma(k)$')

[<matplotlib.lines.Line2D at 0x1fc8ab9cc90>]

Elección de la Best Matching Unit(BMU) y Actualización

Para cada instancia que se presenta al SOM se elige una neurona cuyo vector de referencia cumpla una condición de mínimo, en concreto, de mínima distancia euclídea. Una vez hecha esta elección, los vectores de referencia se van actualizando a medida que las instancias son procesadas por la red auto-organizativa.

Esta actualización se lleva a cabo mediante la fórmula de Kohonen con un núcleo de entorno exponencial.

def BMU(dato,SOM):
    distancias = []
    for ii in np.arange(neu):
        neurona = SOM.loc[ii].to_numpy()
        dist = np.sqrt(np.sum(((dato-neurona)**2)))
        distancias.append(dist)
    # seleccionamos la menor
    
    menorPosicion = np.argmin(np.asarray(distancias))
    return menorPosicion

def actualizarPesos(dato,SOM,posicionBMU,itt):
    '''
    itt: iteracion
    '''
    Vk = []
    for k in np.arange(neu):
        kernel = np.exp(-((posicionBMU-k)**2)/(2*sigma(i)**2))
        SOM.loc[k] = SOM.loc[k]+alpha(itt)*kernel*(dato-SOM.loc[k])
        Vk.append(np.sqrt(np.sum((SOM.loc[posicionBMU]-SOM.loc[k])**2)))
    Vmeanitt = np.mean(Vk)
    return SOM,Vmeanitt

def crearSOM():
    W = np.random.uniform(low=0.0,high=1.0,size=(neu,nd))
    Wdf = pd.DataFrame(W)
    return Wdf

#neu = 25
#ni = 1000
#nd = 10

Vm = []
Tm = []
Tvar = []
Bm = []
Bvar = []
closesdist = []

SOM = crearSOM()
SOMinit = SOM.copy()
for ii in np.arange(ni):
    if ii%100==0:
        print('Etapa:',ii)
    dato = datosIniciales.loc[ii].to_numpy()
    posicionBMU = BMU(dato,SOM)
    cd = np.sqrt(np.sum((dato-SOM.loc[posicionBMU])**2))
    closesdist.append(cd)
    #print(ii,posicionBMU)
    SOM,Vmeanitt = actualizarPesos(dato,SOM,posicionBMU,ii)
    Vm.append(Vmeanitt)
    Tm.append(SOM.to_numpy().ravel().mean())
    Tvar.append(SOM.to_numpy().ravel().std())
    Bm.append(SOM.loc[posicionBMU].to_numpy().mean())
    Bvar.append(SOM.loc[posicionBMU].to_numpy().std())
print('Finished')

Show code cell output Hide code cell output

Etapa: 0

Etapa: 100

Etapa: 200

Etapa: 300

Etapa: 400

Etapa: 500

Etapa: 600

Etapa: 700

Etapa: 800

Etapa: 900

Finished

Comparación de barplots

f,ax=plt.subplots(nrows=1,ncols=2,figsize=(30,15))
SOMinit.plot.bar(stacked=True,ax=ax[0],legend=False)
SOM.plot.bar(stacked=True,ax=ax[1],legend=False)

<Axes: >

Comparación de histogramas

f,ax=plt.subplots(nrows=1,ncols=3,figsize=(30,10))
ax[0].hist(datosIniciales.to_numpy().ravel(),range=[0,1],bins=20)
ax[0].set_title('Datos Iniciales')
ax[1].hist(SOMinit.to_numpy().ravel(),range=[0,1],bins=20)
ax[1].set_title('SOM untrained')
ax[2].hist(SOM.to_numpy().ravel(),range=[0,1],bins=20)
ax[2].set_title('SOM trained')
plt.plot()

[]

Histogramas de un feature

feature  = np.random.choice(nd)
f,ax=plt.subplots(nrows=1,ncols=3,figsize=(30,10))
ax[0].hist(datosIniciales.loc[:,feature].to_numpy().ravel(),range=[0,1],bins=20)
ax[0].set_title('Datos Iniciales')
ax[1].hist(SOMinit.loc[:,feature].to_numpy().ravel(),range=[0,1],bins=20)
ax[1].set_title('SOM untrained')
ax[2].hist(SOM.loc[:,feature].to_numpy().ravel(),range=[0,1],bins=20)
ax[2].set_title('SOM trained')
f.suptitle('Feature: {}'.format(feature))
plt.plot()

[]

f,ax=plt.subplots(nrows=2,ncols=2,figsize=(10,10))
x=np.arange(nd)
ax[0,0].errorbar(x,SOMinit.mean(),SOMinit.std(),linestyle='None', marker='^')
ax[0,0].grid()
ax[0,0].set_title('Por atributo')
ax[1,0].errorbar(x,SOM.mean(),SOM.std(),linestyle='None', marker='^')
ax[1,0].grid()
ax[1,0].set_title('Por atributo')

x = np.arange(neu)
ax[0,1].errorbar(x,SOMinit.mean(axis=1),SOMinit.std(axis=1),linestyle='None', marker='^',alpha=0.1)
ax[0,1].grid()
ax[0,1].set_title('Por neurona')
ax[1,1].errorbar(x,SOM.mean(axis=1),SOM.std(axis=1),linestyle='None', marker='^',alpha=0.1)
ax[1,1].grid()
ax[1,1].set_title('Por neurona')

Text(0.5, 1.0, 'Por neurona')

Asociar datos a neuronas

BMUV = []
for ii in range(ni):
    dato = datosIniciales.loc[ii]
    posicionBMU = BMU(dato,SOM)
    BMUV.append(posicionBMU)

f,ax=plt.subplots(nrows=1,ncols=1,figsize=(8,8))
plt.plot(range(ni),BMUV,'.')
plt.xlabel('Dato')
plt.ylabel('Neurona BMU')

Text(0, 0.5, 'Neurona BMU')

Gráficas de evolucion

datosIniciales.to_numpy().ravel().mean(),datosIniciales.to_numpy().ravel().std()

(0.1668, 0.3727972102899913)

xx = np.arange(ni)
f,axes=plt.subplots(nrows=2,ncols=3,figsize=(30,20))
ax=axes.ravel()
ax[0].plot(xx,Vm,'.')
ax[0].set_title('Vm')
ax[1].plot(xx,Tm,'.')
ax[1].set_title('Tm')
ax[2].plot(xx,Tvar,'.')
ax[2].set_title('Tvar')
ax[3].plot(xx,Bm,'.')
ax[3].set_title('Bm')
ax[4].plot(xx,Bvar,'.')
ax[4].set_title('Bvar')
ax[5].plot(xx,closesdist,'.')
ax[5].set_title('closesdist')

Text(0.5, 1.0, 'closesdist')

Ejercicios

• Realiza un entrenamiento SOM con varias épocas de aprendizaje

• Aleatoriza la entrada de los datos al proceso de aprendizaje

• Varía los parámetros de aprendizaje

• Compara los resultados

• Comenta las gráficas