Detección de anomalías con PCA y python

Detección de anomalías con PCA y python

Joaquín Amat Rodrigo
Diciembre, 2020 (actualizado Octubre 2021)

Introducción


La detección de anomalías (outliers) con Análisis de Componentes Principales (PCA) es una estrategia no supervisada para identificar anomalías cuando los datos no están etiquetados, es decir, no se conoce la clasificación real (anomalía - no anomalía) de las observaciones.

Si bien esta estrategia hace uso de PCA, no utiliza directamente su resultado como forma de detectar anomalías, sino que emplea el error de reconstrucción producido al revertir la reducción de dimensionalidad. El error de reconstrucción como estrategia para detectar anomalías se basa en la siguiente idea: los métodos de reducción de dimensionalidad permiten proyectar las observaciones en un espacio de menor dimensión que el espacio original, a la vez que tratan de conservar la mayor información posible. La forma en que consiguen minimizar la pérdida global de información es buscando un nuevo espacio en el que la mayoría de observaciones puedan ser bien representadas.

El método PCA crea una función que mapea la posición que ocupa cada observación en el espacio original con el que ocupa en el nuevo espacio generado. Este mapeo funciona en ambas direcciones, por lo que también se puede ir desde el nuevo espacio al espacio original. Solo aquellas observaciones que hayan sido bien proyectadas podrán volver a la posición que ocupaban en el espacio original con una precisión elevada.

Dado que la búsqueda de ese nuevo espacio ha sido guiada por la mayoría de las observaciones, serán las observaciones más próximas al promedio las que mejor puedan ser proyectadas y en consecuencia mejor reconstruidas. Las observaciones anómalas, por el contrario, serán mal proyectadas y su reconstrucción será peor. Es este error de reconstrucción (elevado al cuadrado) el que puede emplearse para identificar anomalías.

Para conocer en más detalle el funcionamiento del PCA y su aplicación en python consultar Análisis de Componentes Principales con python.

Librerías


Las librerías utilizadas en este documento son:

In [3]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from mat4py import loadmat

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns
style.use('ggplot') or plt.style.use('ggplot')

# Preprocesado y modelado
# ==============================================================================
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from tqdm import tqdm

import h2o
from h2o.estimators import H2OGeneralizedLowRankEstimator

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')

Datos


Los datos empleados en este documento se han obtenido de Outlier Detection DataSets (ODDS), un repositorio con datos comúnmente empleados para comparar la capacidad que tienen diferentes algoritmos a la hora de identificar anomalías (outliers). Shebuti Rayana (2016). ODDS Library [http://odds.cs.stonybrook.edu]. Stony Brook, NY: Stony Brook University, Department of Computer Science.

Todos estos conjuntos de datos están etiquetados, se conoce si las observaciones son o no anomalías (variable y). Aunque los métodos que se describen en el documento son no supervisados, es decir, no hacen uso de la variable respuesta, conocer la verdadera clasificación permite evaluar su capacidad para identificar correctamente las anomalías.

  • Cardiotocogrpahy dataset link:

    • Número de observaciones: 1831
    • Número de variables: 21
    • Número de outliers: 176 (9.6%)
    • y: 1 = outliers, 0 = inliers
    • Observaciones: todas las variables están centradas y escaladas (media 0, sd 1).
    • Referencia: C. C. Aggarwal and S. Sathe, “Theoretical foundations and algorithms for outlier ensembles.” ACM SIGKDD Explorations Newsletter, vol. 17, no. 1, pp. 24–47, 2015. Saket Sathe and Charu C. Aggarwal. LODES: Local Density meets Spectral Outlier Detection. SIAM Conference on Data Mining, 2016.
  • Speech dataset link:

    • Número de observaciones: 3686
    • Número de variables: 400
    • Número de outliers: 61 (1.65%)
    • y: 1 = outliers, 0 = inliers
    • Referencia: Learing Outlier Ensembles: The Best of Both Worlds – Supervised and Unsupervised. Barbora Micenkova, Brian McWilliams, and Ira Assent, KDD ODD2 Workshop, 2014.
  • Shuttle dataset link:

    • Número de observaciones: 49097
    • Número de variables: 9
    • Número de outliers: 3511 (7%)
    • y: 1 = outliers, 0 = inliers
    • Referencia: Abe, Naoki, Bianca Zadrozny, and John Langford. “Outlier detection by active learning.” Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 2006.

Los datos están disponibles en formato matlab (.mat). Para leer su contenido se emplea la función loadmat() del paquete mat4py 0.1.0.

In [4]:
# Lectura de datos
# ==============================================================================
datos = loadmat(filename='cardio.mat')
datos_X = pd.DataFrame(datos['X'])
datos_X.columns = ["col_" + str(i) for i in datos_X.columns]
datos_y = pd.DataFrame(datos['y'])
datos_y = datos_y.to_numpy().flatten()

PCA


Principal Component Analysis (PCA) es un método estadístico que permite simplificar la complejidad de espacios muestrales con múltiples dimensiones a la vez que conserva su información. Supóngase que existe una muestra con $n$ individuos cada uno con $p$ variables ($X_1$, $X_2$, ..., $X_p$), es decir, el espacio tiene $p$ dimensiones. PCA permite encontrar un número de factores subyacentes, $(z < p)$ que explican aproximadamente lo mismo que las $p$ variables originales. Donde antes se necesitaban $p$ valores para caracterizar a cada individuo, ahora bastan $z$ valores. Cada una de estas $z$ nuevas variables recibe el nombre de componente principal.

Cada componente principal ($Z_i$) se obtiene por combinación lineal de las variables originales. Se pueden entender como nuevas variables obtenidas al combinar de una determinada forma las variables originales. La primera componente principal de un grupo de variables ($X_1$, $X_2$, ..., $X_p$) es la combinación lineal normalizada de dichas variables que tiene mayor varianza. Una vez calculada la primera componente ($Z_1$), se calcula la segunda ($Z_2$) repitiendo el mismo proceso, pero añadiendo la condición de que la combinación lineal no puede estar correlacionada con la primera componente. El proceso se repite de forma iterativa hasta calcular todas las posibles componentes (min(n-1, p)) o hasta que se decida detener el proceso.

La principal limitación que tiene el PCA como método de reducción de dimensionalidad es que solo contempla combinaciones lineales de las variables originales, lo que significa que no es capaz de capturar otro tipo de relaciones.

Modelo


La clase sklearn.decomposition.PCA incorpora las principales funcionalidades que se necesitan a la hora de trabajar con modelos PCA. El argumento n_components determina el número de componentes calculados. Si se indica None, se calculan todas las posibles (min(filas, columnas) - 1).

Por defecto, PCA() centra los valores pero no los escala. Esto es importante ya que, si las variables tienen distinta dispersión, como en este caso, es necesario escalarlas. Una forma de hacerlo es combinar un StandardScaler() y un PCA() dentro de un pipeline. Para más información sobre el uso de pipelines consultar Pipeline y ColumnTransformer.

In [5]:
# Entrenamiento modelo PCA con escalado de los datos
# ==============================================================================
pca_pipeline = make_pipeline(StandardScaler(), PCA(n_components=None))
pca_pipeline.fit(X=datos_X)

# Se extrae el modelo entrenado del pipeline para su diagnóstico
modelo_pca = pca_pipeline.named_steps['pca']

Diagnóstico


Tras obtener el PCA, se evalúa su capacidad para reducir la dimensionalidad de los datos explorando la varianza explicada y acumulada de cada componente.

In [6]:
# Porcentaje de varianza explicada por cada componente
# ==============================================================================
print('----------------------------------------------------')
print('Porcentaje de varianza explicada por cada componente')
print('----------------------------------------------------')
print(modelo_pca.explained_variance_ratio_)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 4))
ax.bar(
    x      = np.arange(modelo_pca.n_components_) + 1,
    height = modelo_pca.explained_variance_ratio_
)

for x, y in zip(np.arange(modelo_pca.n_components_) + 1, modelo_pca.explained_variance_ratio_):
    label = round(y, 2)
    ax.annotate(
        label,
        (x,y),
        textcoords="offset points",
        xytext=(0,10),
        ha='center'
    )

ax.set_xticks(np.arange(modelo_pca.n_components_) + 1)
ax.set_ylim(0, 1.1)
ax.set_title('Porcentaje de varianza explicada por cada componente')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza explicada');
----------------------------------------------------
Porcentaje de varianza explicada por cada componente
----------------------------------------------------
[2.69831338e-01 1.75787255e-01 8.68938668e-02 6.88155848e-02
 5.95308477e-02 4.97905270e-02 4.58906579e-02 4.42857769e-02
 3.75980654e-02 3.19494743e-02 2.92016595e-02 2.56986810e-02
 1.92184240e-02 1.74753805e-02 1.31181188e-02 8.86498538e-03
 6.71579614e-03 5.51775331e-03 2.43704462e-03 1.37876279e-03
 5.19122572e-33]
In [7]:
# Porcentaje de varianza explicada acumulada
# ==============================================================================
prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
print('------------------------------------------')
print('Porcentaje de varianza explicada acumulada')
print('------------------------------------------')
print(prop_varianza_acum)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 4))
ax.plot(
    np.arange(modelo_pca.n_components_) + 1,
    prop_varianza_acum,
    marker = 'o'
)

for x, y in zip(np.arange(modelo_pca.n_components_) + 1, prop_varianza_acum):
    label = round(y, 2)
    ax.annotate(
        label,
        (x,y),
        textcoords="offset points",
        xytext=(0,10),
        ha='center'
    )

ax.axvline(x=11, linestyle = '--')
ax.set_ylim(0, 1.1)
ax.set_xticks(np.arange(modelo_pca.n_components_) + 1)
ax.set_title('Porcentaje de varianza explicada acumulada')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza acumulada');
------------------------------------------
Porcentaje de varianza explicada acumulada
------------------------------------------
[0.26983134 0.44561859 0.53251246 0.60132804 0.66085889 0.71064942
 0.75654008 0.80082585 0.83842392 0.87037339 0.89957505 0.92527373
 0.94449216 0.96196754 0.97508566 0.98395064 0.99066644 0.99618419
 0.99862124 1.         1.        ]

Con las primeras 11 componentes se consigue explicar aproximadamente el 90% de la varianza observada en los datos.

Reconstrucción


Una vez obtenido el PCA (matriz de componentes (eigenvectors), proyecciones y medias), la reconstrucción de las observaciones iniciales se puede obtener empleando la siguiente ecuación:

$$\text{reconstrucción} = \text{proyecciones} \cdot \text{componentes}^\top$$

Es importante tener en cuenta que, si los datos han sido centrados y escalados (cosa que en términos generales debe hacerse al aplicar PCA) esta transformación debe revertirse.

Con el método inverse_transform() de la clase PCA puede revertirse la transformación y reconstruir el valor inicial. Es importante tener en cuenta que, la reconstrucción, solo será completa si se han incluido todas las componentes. El método inverse_transform() utiliza para la reconstrucción todas las componentes incluidas al crear el objeto PCA.

In [8]:
# Reconstrucción de las proyecciones
# ==============================================================================
proyecciones = pca_pipeline.transform(datos_X)
reconstruccion = pca_pipeline.inverse_transform(X=proyecciones)
reconstruccion = pd.DataFrame(
                    reconstruccion,
                    columns = datos_X.columns,
                    index   = datos_X.index
)
print('------------------')
print('Valores originales')
print('------------------')
display(reconstruccion.head(3))

print('')
print('---------------------')
print('Valores reconstruidos')
print('---------------------')
display(datos_X.head(3))
------------------
Valores originales
------------------
col_0 col_1 col_2 col_3 col_4 col_5 col_6 col_7 col_8 col_9 ... col_11 col_12 col_13 col_14 col_15 col_16 col_17 col_18 col_19 col_20
0 0.004912 0.693191 -0.20364 0.595322 0.353190 -0.061401 -0.278295 -1.650444 0.759072 -0.420487 ... 1.485973 -0.798376 1.854728 0.622631 0.963083 0.301464 0.193113 0.231498 -0.289786 -0.493294
1 0.110729 -0.079903 -0.20364 1.268942 0.396246 -0.061401 -0.278295 -1.710270 0.759072 -0.420487 ... 1.485973 -0.798376 1.854728 0.278625 0.963083 0.301464 0.129265 0.093563 -0.256385 -0.493294
2 0.216546 -0.272445 -0.20364 1.050988 0.148753 -0.061401 -0.278295 -1.710270 1.106509 -0.420487 ... 1.141780 -1.332931 0.314688 2.342663 -0.488279 0.061002 0.065417 0.024596 -0.256385 1.140018

3 rows × 21 columns

---------------------
Valores reconstruidos
---------------------
col_0 col_1 col_2 col_3 col_4 col_5 col_6 col_7 col_8 col_9 ... col_11 col_12 col_13 col_14 col_15 col_16 col_17 col_18 col_19 col_20
0 0.004912 0.693191 -0.20364 0.595322 0.353190 -0.061401 -0.278295 -1.650444 0.759072 -0.420487 ... 1.485973 -0.798376 1.854728 0.622631 0.963083 0.301464 0.193113 0.231498 -0.289786 -0.493294
1 0.110729 -0.079903 -0.20364 1.268942 0.396246 -0.061401 -0.278295 -1.710270 0.759072 -0.420487 ... 1.485973 -0.798376 1.854728 0.278625 0.963083 0.301464 0.129265 0.093563 -0.256385 -0.493294
2 0.216546 -0.272445 -0.20364 1.050988 0.148753 -0.061401 -0.278295 -1.710270 1.106509 -0.420487 ... 1.141780 -1.332931 0.314688 2.342663 -0.488279 0.061002 0.065417 0.024596 -0.256385 1.140018

3 rows × 21 columns

Se comprueba que, utilizando todas las componentes, la reconstrucción es total. Los valores reconstruidos son iguales a los datos originales (las pequeñas diferencias se deben a la imprecisión numérica de las operaciones).

Error de reconstrucción


El error cuadrático medio de reconstrucción de una observación se calcula como el promedio de las diferencias al cuadrado entre el valor original de sus variables y el valor reconstruido, es decir, el promedio de los errores de reconstrucción de todas sus variables elevados al cuadrado.

En el apartado anterior, se mostró que, empleando todas las componentes, la reconstrucción es total, por lo que el error de reconstrucción es cero. Cuando el objetivo de la reconstrucción es identificar observaciones anómalas, hay que emplear únicamente un subconjunto de todas las componentes, permitiendo así diferenciar las observaciones fáciles de reconstruir de las difíciles.

Para facilitar el proceso, se define una función que, dado un conjunto de datos y un número de componentes: escale los datos, entrene un modelo PCA y devuelva la reconstrucción de los datos iniciales y su error.

In [9]:
def reconstruccion_pca(X, n_components, X_new=None):
    '''
    Función para calcular la reconstrucción y error de un conjunto de datos
    empleando un PCA
    
    Parameters
    ----------
    
    X (data.frame): datos de entrenamiento del PCA.
    
    X_new (data.frame): datos sobre los que aplicar la reconstrucción. Si es None,
                        se emplea X.
    
    id_components (list, numpy.narray): indice de las componentes empleadas para
                                        la reconstrucción, empezando por 0.
                                        Por defecto se emplean todas.
    
    Returns
    -------
    
    reconstruccion (data.frame): reconstrucción de los datos.
    
    error_reconstruccion (numpy.narray): error de reconstrucción de las observaciones.
    
    '''
    
    if X_new is None:
        X_new = X
    
    # Entrenamiento modelo PCA con escalado de los datos
    pca_pipeline = make_pipeline(StandardScaler(), PCA(n_components=n_components))
    pca_pipeline.fit(X=X)
    
    # Proyectar los datos
    proyecciones = pca_pipeline.transform(X_new)
    
    # Reconstrucción
    reconstruccion = pca_pipeline.inverse_transform(X=proyecciones)
    reconstruccion = pd.DataFrame(
                        reconstruccion,
                        columns = X_new.columns,
                        index   = X_new.index
                    )
    
    # Error cuadrático medio de reconstrucción
    error_reconstruccion = reconstruccion - X_new
    error_reconstruccion = error_reconstruccion**2
    error_reconstruccion = error_reconstruccion.mean(axis=1)
    
    print(f"Reconstrucción con un PCA de {n_components} componentes")
    
    return reconstruccion, error_reconstruccion
In [10]:
# Reconstrucción con las 11 primeras componentes (90% de la varianza explicada)
# ==============================================================================
reconstruccion, error_reconstruccion = reconstruccion_pca(X=datos_X, n_components=11)
Reconstrucción con un PCA de 11 componentes
In [11]:
# Distribución del error de reconstrucción
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 3.5))
sns.kdeplot(
    error_reconstruccion,
    fill    = True,
    ax      = ax
)
sns.rugplot(error_reconstruccion,  ax=ax, color='black')
ax.set_title('Distribución de los errores de reconstrucción (PCA)')
ax.set_xlabel('Error de reconstrucción');

Detección de anomalías


Una vez que el error de reconstrucción ha sido calculado, se puede emplear como criterio para identificar anomalías. Asumiendo que la reducción de dimensionalidad se ha realizado de forma que la mayoría de los datos (los normales) queden bien representados, aquellas observaciones con mayor error de reconstrucción deberían ser las más atípicas.

En la práctica, si se está empleando esta estrategia de detección es porque no se dispone de datos etiquetados, es decir, no se conoce qué observaciones son realmente anomalías. Sin embargo, como en este ejemplo se dispone de la clasificación real, se puede verificar si realmente los datos anómalos tienen errores de reconstrucción más elevados.

In [12]:
df_resultados = pd.DataFrame({
                    'error_reconstruccion' : error_reconstruccion,
                    'anomalia'             : datos_y
                })

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 3.5))
sns.boxplot(
    x     = 'anomalia',
    y     = 'error_reconstruccion',
    data  = df_resultados,
    #color = "white",
    palette = 'tab10',
    ax    = ax
)
ax.set_yscale("log")
ax.set_title('Distribución de los errores de reconstrucción (PCA)')
ax.set_ylabel('log(Error de reconstrucción)')
ax.set_xlabel('clasificación (0 = normal, 1 = anomalía)');

La distribución de los errores de reconstrucción en el grupo de las anomalías (1) es claramente superior. Sin embargo, al existir solapamiento, si se clasifican las n observaciones con mayor error de reconstrucción como anomalías, se incurriría en errores de falsos positivos.

Acorde a la documentación, el set de datos Cardiotocogrpahy contiene 176 anomalías. Véase la matriz de confusión resultante si se clasifican como anomalías las 176 observaciones con mayor error de reconstrucción.

In [13]:
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = df_resultados \
                .sort_values('error_reconstruccion', ascending=False) \
                .reset_index(drop=True)

df_resultados['clasificacion'] = np.where(df_resultados.index <= 176, 1, 0)

pd.crosstab(
    df_resultados.anomalia,
    df_resultados.clasificacion
)
Out[13]:
clasificacion 0 1
anomalia
0.0 1545 110
1.0 109 67

De las 176 observaciones identificadas como anomalías, solo el 38% (67/176) lo son. El porcentaje de falsos positivos (62%) es muy elevado.

Últimas componentes


En el apartado anterior, se han empleado las 11 primeras componentes de un total de 21 ya que con ellas se puede explicar aproximadamente el 90% de la varianza observada. Que esta selección de componentes sea adecuada para la detección de anomalías depende en gran medida de en qué direcciones se desvíen las anomalías. Con frecuencia, los datos anómalos no tienen valores atípicos en las direcciones dominantes pero sí en las poco dominantes, por lo que su comportamiento atípico quedaría recogido en las últimas componentes. Se repite la detección de anomalías pero esta vez descartando las 4 primeras componentes principales en la reconstrucción.

El método inverse_transform() de la clase PCA no permite seleccionar componentes específicas con las que reconstruir los datos. Afortunadamente, la reconstrucción puede hacerse fácilmente multiplicado las proyecciones por los vectores que definen cada componente.

In [14]:
def reconstruccion_pca(X, X_new=None, id_components=None):
    '''
    Función para calcular la reconstrucción y error de un conjunto de datos
    empleando un un subconjunto de componentes PCA
    
    Parameters
    ----------
    
    X (data.frame): datos de entrenamiento del PCA.
    
    X_new (data.frame): datos sobre los que aplicar la reconstrucción. Si es None,
                        se emplea X.
    
    id_components (list, numpy.narray): indice de las componentes empleadas para
                                        la reconstrucción, empezando por 0.
                                        Por defecto se emplean todas.
    
    Returns
    -------
    
    reconstruccion (data.frame): reconstrucción de los datos.
    
    error_reconstruccion (numpy.narray): error de reconstrucción de cada observación.
    '''
    
    if X_new is None:
        X_new = X
    
    # Si se especifica id_components, se calculan solo tantas como necesarias.
    if id_components is None:
        n_components = None
        id_components = np.arange(min(X.shape))
    else:
        n_components = max(id_components)
        id_components = np.array(id_components)
    
    # Entrenamiento modelo PCA con escalado de los datos
    pca_pipeline = make_pipeline(StandardScaler(), PCA(n_components=n_components))
    pca_pipeline.fit(X=X)
    
    # Proyectar los datos
    proyecciones = pca_pipeline.transform(X_new)
    
    # Reconstrucción
    components = pca_pipeline['pca'].components_
    
    reconstruccion = np.dot(
                        proyecciones[:, id_components-1],
                        components[id_components-1, :]
                     )
    
    # Se deshace el escalado
    reconstruccion = pca_pipeline['standardscaler'].inverse_transform(reconstruccion)
    
    reconstruccion = pd.DataFrame(
                        reconstruccion,
                        columns = X_new.columns,
                        index   = X_new.index
                    )
    
    # Error cuadrático medio de reconstrucción
    error_reconstruccion = reconstruccion - X_new
    error_reconstruccion = error_reconstruccion**2
    error_reconstruccion = error_reconstruccion.mean(axis=1)
    
    print(f"Reconstrucción con las componentes: {id_components}")
    
    return reconstruccion, error_reconstruccion
In [15]:
# Reconstrucción empleando todas excepto las 4 primeras componentes
# ==============================================================================
reconstruccion, error_reconstruccion = reconstruccion_pca(X=datos_X, id_components = range(4,22))
Reconstrucción con las componentes: [ 4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21]
In [16]:
df_resultados = pd.DataFrame({
                    'error_reconstruccion' : error_reconstruccion,
                    'anomalia'             : datos_y
                })

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 3.5))
sns.boxplot(
    x     = 'anomalia',
    y     = 'error_reconstruccion',
    data  = df_resultados,
    #color = "white",
    palette = 'tab10',
    ax    = ax
)
ax.set_yscale("log")
ax.set_title('Distribución de los errores de reconstrucción (PCA)')
ax.set_ylabel('log(Error de reconstrucción)')
ax.set_xlabel('clasificación (0 = normal, 1 = anomalía)');
In [17]:
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = df_resultados \
                .sort_values('error_reconstruccion', ascending=False) \
                .reset_index(drop=True)

df_resultados['clasificacion'] = np.where(df_resultados.index <= 176, 1, 0)

pd.crosstab(
    df_resultados.anomalia,
    df_resultados.clasificacion
)
Out[17]:
clasificacion 0 1
anomalia
0.0 1601 54
1.0 53 123

Descartando las primeras 4 componentes, 123 de las 176 (70%%) observaciones identificadas como anomalías lo son realmente. Se ha conseguido reducir el porcentaje de falsos positivos al 30%.

Reentrenamiento iterativo


El PCA anterior se ha obtenido empleando todas las observaciones, incluyendo potenciales anomalías. Dado que el objetivo es generar un espacio de proyección únicamente con los datos “normales”, se puede mejorar el resultado reentrenando el PCA pero excluyendo las n observaciones con mayor error de reconstrucción (potenciales anomalías).

Se repite la detección de anomalías pero, esta vez, descartando las observaciones con un error de reconstrucción superior al cuantil 0.9.

In [18]:
# Se eliminan las observaciones con un error de reconstrucción superior al
# cuantil 0.9
# ==============================================================================
cuantil = np.quantile(a=error_reconstruccion,  q=0.9)
datos_X_trimmed = datos_X.loc[error_reconstruccion < cuantil, :].copy()
datos_y_trimmed = datos_y[error_reconstruccion < cuantil].copy()
In [19]:
# Reconstrucción empleando todas excepto las 4 primeras componentes
# ==============================================================================
reconstruccion2, error_reconstruccion2 = reconstruccion_pca(
                                            X     = datos_X_trimmed,
                                            X_new = datos_X,
                                            id_components = range(4,21)
                                        )
Reconstrucción con las componentes: [ 4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20]
In [20]:
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = pd.DataFrame({
                    'error_reconstruccion' : error_reconstruccion2,
                    'anomalia'             : datos_y
                })

df_resultados = df_resultados \
                .sort_values('error_reconstruccion', ascending=False) \
                .reset_index(drop=True)

df_resultados['clasificacion'] = np.where(df_resultados.index <= 176, 1, 0)

pd.crosstab(
    df_resultados.anomalia,
    df_resultados.clasificacion
)
Out[20]:
clasificacion 0 1
anomalia
0.0 1605 50
1.0 49 127

Tras descartar el 10% de las observaciones con mayor error y reentrenando el PCA, 127 de las 176 (72%) observaciones identificadas como anomalías lo son realmente. Se ha reducido ligeramente el porcentaje de falsos positivos al 28%.

Bagging PCA


El método de bagging (bootstrap aggregating) permite reducir la varianza de modelos estadísticos y de machine learning. En términos generales, este método consiste entrenar n modelos, cada uno con una pseudo-muestra obtenida por resampling a partir de los datos de entrenamiento. La predicción final de una observación es el promedio de las predicciones de todos los modelos generados.

Esta estrategia puede combinarse con PCA para la detección de anomalías. En este caso, se entrenan múltiples PCA, cada uno con un resampling de los datos iniciales. Con cada uno de los PCA se calcula el error de reconstrucción de las observaciones originales. Finalmente, se promedia el error de reconstrucción de cada muestra.

Funciones

In [21]:
# Funciones
# ==============================================================================

def reconstruccion_pca(X, X_new=None, id_components=None):
    '''
    Función para entrenar un modelo PCA y obtener a partir de este la reconstrucción
    de nuevas observaciones y su error asociado.
    
    Parameters
    ----------
    
    X : pd.DataFrame, np.narray
        Datos de entrenamiento del PCA.
    
    X_new : pd.DataFrame, np.narray, default `None`
        Datos sobre los que aplicar la reconstrucción. Si es `None`, se emplea `X`.
    
    id_components : list, np.narray, default `None`
        Índice de las componentes empleadas para la reconstrucción, empezando por 0.
        Si es `None`, se emplean todas.
    
    Returns
    -------
    
    reconstruccion: pd.DataFrame
        Reconstrucción de los datos.
    
    error_reconstruccion : np.narray
        Error de reconstrucción de cada observación.
    '''
    
    if X_new is None:
        X_new = X
    
    # Si se especifica id_components, se calculan solo tantas como necesarias.
    if id_components is None:
        n_components = None
        id_components = np.arange(min(X.shape))
    else:
        n_components = max(id_components)
        id_components = np.array(id_components)
    
    # Entrenamiento modelo PCA con escalado de los datos
    pca_pipeline = make_pipeline(StandardScaler(), PCA(n_components=n_components))
    pca_pipeline.fit(X=X)
    
    # Proyectar los datos
    proyecciones = pca_pipeline.transform(X_new)
    
    # Reconstrucción
    components = pca_pipeline['pca'].components_
    
    reconstruccion = np.dot(
                        proyecciones[:, id_components-1],
                        components[id_components-1, :]
                     )
    
    # Se deshace el escalado
    reconstruccion = pca_pipeline['standardscaler'].inverse_transform(reconstruccion)
    
    reconstruccion = pd.DataFrame(
                        reconstruccion,
                        columns = X_new.columns,
                        index   = X_new.index
                    )
    
    # Error cuadrático medio de reconstrucción
    error_reconstruccion = reconstruccion - X_new
    error_reconstruccion = error_reconstruccion**2
    error_reconstruccion = error_reconstruccion.mean(axis=1)
        
    return reconstruccion, error_reconstruccion


def bagging_reconstruccion_pca(X, X_new=None, iteraciones=500, verbose=0):
    '''
    Función para calcular el error de reconstrucción combinando el método de
    bagging con el de PCA. En cada iteración, se entrena un modelo PCA con
    un resample de los datos y se reconstruye con un conjunto aleatorio de
    componentes. El conjunto de componentes puede tener entre 1 y la mitad de
    las componentes máximas.
    
    Parameters
    ----------
    
    X (data.frame): datos de entrenamiento del PCA.
    
    X_new (data.frame): datos sobre los que aplicar la reconstrucción. Si es None,
                        se emplea X.
    
    iteraciones (int): número de iteraciones de bagging.
    
    Returns
    -------
    
    error_reconstruccion (numpy.narray): error de reconstrucción promedio de
                                         cada observación.
    '''
    
    if X_new is None:
        X_new = X
        
    mat_error = np.full(shape=(X_new.shape[0], iteraciones), fill_value=np.nan)
    
    for i in tqdm(range(iteraciones)):
        
        # Resampling
        # ----------------------------------------------------------------------        
        resample = X.sample(frac=0.632, axis=0)
        resample = X.sample(frac=0.8, axis=1)
        
        # Si debido al resampling, alguna columna es constante, se elimina
        if any(resample.var() == 0):
            constat_cols = resample.columns[resample.var() == 0]
            resample = resample.drop(columns=constat_cols)
            
            
        # Selección aleatoria de las componentes empleadas en la reconstrucción
        # ----------------------------------------------------------------------
        n_componentes = min(resample.shape)
        id_componentes = np.random.choice(
                            range(n_componentes),
                            size = max(1, int(n_componentes/2)),
                            replace = False
                         )
        
        if verbose > 0:
            print(f"Componentes utilizadas iteración {i}: {np.sort(id_componentes)}")
            
        # Error reconstrucción
        # ----------------------------------------------------------------------            
        _, error_reconstruccion = reconstruccion_pca(
                                    X = resample,
                                    X_new = X_new.loc[:, resample.columns],
                                    id_components=id_componentes
                                  )   
        
        # Almacenamiento resultados iteración
        # ----------------------------------------------------------------------
        mat_error[:, i] = error_reconstruccion
        
    # Valor promedio de todas las iteraciones (se ecalan los valores antes)
    # ----------------------------------------------------------------------
    mat_error = StandardScaler().fit_transform(mat_error)
    error_promedio = np.nanmean(mat_error, axis=1)

    return error_promedio

Error de reconstrucción

In [22]:
# Cálculo error de reconstrucción
# ==============================================================================
error_reconstruccion = bagging_reconstruccion_pca(X=datos_X, iteraciones=500)
100%|█████████████████████████████████████████| 500/500 [00:15<00:00, 32.71it/s]

Detección de anomalías

In [23]:
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = pd.DataFrame({
                    'error_reconstruccion' : error_reconstruccion,
                    'anomalia'             : datos_y
                })

df_resultados = df_resultados \
                .sort_values('error_reconstruccion', ascending=False) \
                .reset_index(drop=True)

df_resultados['clasificacion'] = np.where(df_resultados.index <= 176, 1, 0)

pd.crosstab(
    df_resultados.anomalia,
    df_resultados.clasificacion
)
Out[23]:
clasificacion 0 1
anomalia
0.0 1582 73
1.0 72 104

La combinación de Bagging y PCA consigue unos resultados notablemente superiores a los que se consigue únicamente con PCA (sin filtrado iterativo de observaciones).

Generalized Low Rank Model (GLRM)

Una de las limitaciones del método de PCA para su uso en la detección de anomalías es su alta sensibilidad a la presencia de valores anómalos (outliers). Dado que el objetivo es aprender un nuevo espacio en el que la mayoría de observaciones (las normales) queden bien representadas, conviene minimizar la influencia de los outliers presentes en el conjunto de entrenamiento.

Los modelos generalized low rank models (GLRMs) son una generalización del PCA y de la factorización matricial que, entre otras características, permite utilizar funciones de coste más robustas que el error cuadrático (empleado por el PCA), por ejemplo, el error absoluto o Huber. Son, por lo tanto, un método de reducción de dimensionalidad potencialmente más adecuado que el PCA para ser utilizado en la detección de anomalías.

Desde python, es posible acceder a este tipo de modelos gracias a la implementación disponible en la librería H2O. Para más información sobre cómo utilizar esta excelente librería consultar Machine learning con H2O y Python.

Nota: un GLRM con el error absoluto como función de coste es equivalente a un PCA-Robusto.

Modelo

In [33]:
# Iniciación del cluster y tranferencia de datos
# ==============================================================================
h2o.init()
h2o.remove_all()
h2o.no_progress()

datos_X_h2o = h2o.H2OFrame(datos_X)
Checking whether there is an H2O instance running at http://localhost:54321 . connected.
H2O_cluster_uptime: 31 secs
H2O_cluster_timezone: Europe/Madrid
H2O_data_parsing_timezone: UTC
H2O_cluster_version: 3.34.0.1
H2O_cluster_version_age: 17 days
H2O_cluster_name: H2O_from_python_ximo_xo3xv7
H2O_cluster_total_nodes: 1
H2O_cluster_free_memory: 1.879 Gb
H2O_cluster_total_cores: 8
H2O_cluster_allowed_cores: 8
H2O_cluster_status: locked, healthy
H2O_connection_url: http://localhost:54321
H2O_connection_proxy: {"http": null, "https": null}
H2O_internal_security: False
H2O_API_Extensions: Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4
Python_version: 3.7.9 final
In [25]:
# Entrenamiento modelo GLRM con el máximo de componenetes posibles
# ==============================================================================
glrm_model = H2OGeneralizedLowRankEstimator(
                k         = min(datos_X_h2o.shape),
                loss      = "absolute",
                transform = 'standardize'
            )
_ = glrm_model.train(training_frame=datos_X_h2o)
In [26]:
# Porcentaje de varianza acumulada
# ==============================================================================
prop_varianza_acum = glrm_model.varimp(use_pandas=True).iloc[2, 1:]
print('------------------------------------------')
print('Porcentaje de varianza explicada acumulada')
print('------------------------------------------')
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9, 4))
ax.plot(
    np.arange(len(prop_varianza_acum)) + 1,
    prop_varianza_acum,
    marker = 'o'
)

for x, y in zip(np.arange(len(prop_varianza_acum)) + 1, prop_varianza_acum):
    label = round(y, 2)
    ax.annotate(
        label,
        (x,y),
        textcoords="offset points",
        xytext=(0,10),
        ha='center'
    )

ax.set_ylim(0, 1.1)
ax.set_xticks(np.arange(modelo_pca.n_components_) + 1)
ax.set_title('Porcentaje de varianza explicada acumulada')
ax.set_xlabel('Componente')
ax.set_ylabel('Por. varianza acumulada');
------------------------------------------
Porcentaje de varianza explicada acumulada
------------------------------------------
In [27]:
glrm_model = H2OGeneralizedLowRankEstimator(
                k         = 14,
                loss      = "absolute",
                transform = 'standardize'
            )
_ = glrm_model.train(training_frame=datos_X_h2o)

Proyección y reconstrucción


Con el método reconstruct de un modelo H2OGeneralizedLowRankEstimator se puede obtener directamente la reconstrucción para cada observación. Internamente, este método está obteniendo primero las proyecciones para cada observación y luego revirtiendo el proceso.

In [28]:
reconstruccion = glrm_model.reconstruct(test_data=datos_X_h2o, reverse_transform=True)

Una vez obtenidas las reconstrucciones, se puede calcular el error cuadrático medio de reconstrucción de una observación con el promedio de las diferencias al cuadrado entre el valor original de sus variables y el valor reconstruido, es decir, el promedio de los errores de reconstrucción de todas sus variables elevados al cuadrado.

In [29]:
error_reconstruccion = reconstruccion - datos_X_h2o
error_reconstruccion = error_reconstruccion.as_data_frame()
error_reconstruccion = error_reconstruccion**2
error_reconstruccion = error_reconstruccion.mean(axis=1)

Detección de anomalías

In [30]:
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = pd.DataFrame({
                    'error_reconstruccion' : error_reconstruccion,
                    'anomalia'             : datos_y
                })

df_resultados = df_resultados \
                .sort_values('error_reconstruccion', ascending=False) \
                .reset_index(drop=True)

df_resultados['clasificacion'] = np.where(df_resultados.index <= 176, 1, 0)

pd.crosstab(
    df_resultados.anomalia,
    df_resultados.clasificacion
)
Out[30]:
clasificacion 0 1
anomalia
0.0 1528 127
1.0 126 50

En contra a lo que cabría esperar, en este caso, un GLRM con función de coste absoluta no consigue mejores resultados que un PCA.

Información de sesión

In [31]:
import session_info
session_info.show(html=False)
-----
h2o                 3.34.0.1
mat4py              0.5.0
matplotlib          3.4.3
numpy               1.19.5
pandas              1.2.5
seaborn             0.11.0
session_info        1.0.0
sklearn             0.24.2
tqdm                4.62.1
-----
IPython             7.26.0
jupyter_client      6.1.7
jupyter_core        4.6.3
jupyterlab          2.1.3
notebook            6.4.0
-----
Python 3.7.9 (default, Aug 31 2020, 12:42:55) [GCC 7.3.0]
Linux-5.11.0-37-generic-x86_64-with-debian-bullseye-sid
-----
Session information updated at 2021-10-02 15:22

Bibliografía


Outlier Analysis Aggarwal, Charu C.

Outlier Ensembles: An Introduction by Charu C. Aggarwal, Saket Sathe

Introduction to Machine Learning with Python: A Guide for Data Scientists

Python Data Science Handbook by Jake VanderPlas

¿Cómo citar este documento?

Detección de anomalías con PCA y python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py21-deteccion-anomalias-pca-python.html DOI


¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


Creative Commons Licence
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.