Detección de anomalías con PCA y python
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.
# 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.
# 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.
# 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.
# 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 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');
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.
# 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))
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.
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
# Reconstrucción con las 11 primeras componentes (90% de la varianza explicada)
# ==============================================================================
reconstruccion, error_reconstruccion = reconstruccion_pca(X=datos_X, n_components=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.
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.
# 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
)
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.
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
# Reconstrucción empleando todas excepto las 4 primeras componentes
# ==============================================================================
reconstruccion, error_reconstruccion = reconstruccion_pca(X=datos_X, id_components = range(4,22))
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)');
# 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
)
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.
# 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()
# 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)
)
# 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
)
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¶
# 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¶
# Cálculo error de reconstrucción
# ==============================================================================
error_reconstruccion = bagging_reconstruccion_pca(X=datos_X, iteraciones=500)
Detección de anomalías¶
# 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
)
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¶
# Iniciación del cluster y tranferencia de datos
# ==============================================================================
h2o.init()
h2o.remove_all()
h2o.no_progress()
datos_X_h2o = h2o.H2OFrame(datos_X)
# 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)
# 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');
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.
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.
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¶
# 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
)
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¶
import session_info
session_info.show(html=False)
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 
¿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! 😊
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.