Detección de anomalías con Gaussian Mixture Model (GMM) y python
Introducción¶
Gaussian Mixture Model (GMM) es un modelo probabilístico en el que se considera que las observaciones siguen una distribución probabilística formada por la combinación de múltiples distribuciones normales (componentes). Puede entenderse como una generalización de K-means con la que, en lugar de asignar cada observación a un único cluster, se obtiene una distribución probabilidad de pertenencia a cada uno.
Ajustar un modelo GMM consiste en estimar los parámetros que definen la función de distribución de cada componente: la media y la matriz de covarianza. Por ejemplo, si el modelo tiene dos componentes, hay que encontrar la media y la matriz de covarianzas de cada una. Si es un problema multidimensional, por ejemplo de 3 variables, la media de cada componente es un vector de 3 valores y la matriz de covarianza una matriz de 3x3. El algoritmo más empleado para realizar el ajuste es Expectation-Maximization (EM).
Una vez aprendidos los parámetros, se puede calcular la densidad de probabilidad que tiene cada observación de pertenecer a cada componente y al conjunto de la distribución. Observaciones con muy poca densidad de probabilidad pueden considerarse como anomalías.
Densidad de probabilidad
Al tratarse de un modelo probabilístico, el ajuste de un modelo GMM genera es en realidad una función de densidad de probabilidad. El concepto de densidad puede entenderse como un análogo al de probabilidad en distribuciones discretas, pero, en lugar de acotada en el rango [0,1], puede tomar valores [0, +$\inf$]. El valor de densidad de probabilidad es una medida relativa de verosimulitud (likelihood), cuanto mayor es el valor de densidad de una observación, mayor evidenca de que la observación pertenece a una determinada distribución.
Con frecuencia, para facilitar los cálculos, en lugar de utilizar el valor de densidad se utiliza el su logaritmo. La interpretación es la misma, cuanto mayor su valor, mayor la evidencia de que la observación pertenece a la distribución.
Aspectos prácticos
Al entrenar un modelo GMM, junto con el número de componentes, hay que determinar el tipo de matriz de covarianza. Dependiendo esto, la forma de las distribuciones de las componentes puede ser:
tied: todas las componentes comparten la misma matriz de covarianza.
diagonal: las dimensiones de cada componente a lo largo de cada dimensión puede ser distintas, pero siempre quedan alineadas con los ejes, es decir, su orientaciones son limitadas.
spherical: las dimensiones de cada componente son las mismas en todas las dimensiones. Esto permite generar distribuciones de distinto tamaño pero todas esféricas.
full: cada componente tiene su propia matriz de covarianza, por lo que pueden tener cualquier orientación y dimensiones.
Gaussian Mixture Model (GMM) en Python
Una de las principales implementaciones de Gaussian Mixture Model está disponibles en Scikit Learn
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from mat4py import loadmat
from sklearn.datasets import make_blobs
# 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.mixture import GaussianMixture
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
# 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.
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()
Ejemplo 1¶
Modelo¶
Con la clase sklearn.mixture.GaussianMixture de Scikit-Learn se pueden entrenar modelos GMMs utilizando el algoritmo expectation-maximization (EM) . Entre sus parámetros destacan:
n_components: número de componentes (en este caso clusters) que forman el modelo.covariance_type: tipo de matriz de covarianza (‘full’ (default), ‘tied’, ‘diag’, ‘spherical’).max_iter: número máximo de iteraciones permitidas.random_state: semilla para garantizar la reproducibilidad de los resultados.
# Definición y entrenamiento del modelo GMM
# ==============================================================================
modelo_gmm = GaussianMixture(
n_components = 4,
covariance_type = 'full',
random_state = 123
)
modelo_gmm.fit(X=datos_X)
El objeto devuelto por GaussianMixture contiene entre otros datos: el peso de cada componente (cluster) en el modelo (weights_), su media (means_) y matriz de covarianza (covariances_). La estructura de esta última depende del tipo de matriz empleada en el ajuste del modelo (covariance_type).
¿Cómo saber el número de componentes y tipo de matriz de covarianza?
Al tratarse de un problema no supervisado, no hay forma de conocer de antemano el número de componentes y tipo de matriz de covarianza óptimos. Afortunadamente, al ser un modelo probabilístico, se puede recurrir a métricas como el Akaike information criterion (AIC) o Bayesian information criterion (BIC) para identificar cómo de bien se ajustan los datos observados al modelo creado, a la vez que se controla el exceso de overfitting. En la implementación de Scikit Learn, para ambas métricas, cuanto más bajo el valor, mejor.
# Tunning del modelo GMM
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
n_components = range(1, 40)
covariance_types = ['spherical', 'tied', 'diag', 'full']
for covariance_type in covariance_types:
valores_bic = []
for i in n_components:
modelo = GaussianMixture(n_components=i, covariance_type=covariance_type)
modelo = modelo.fit(datos_X)
valores_bic.append(modelo.bic(datos_X))
ax.plot(n_components, valores_bic, label=covariance_type)
ax.set_title("Valores BIC")
ax.set_xlabel("Número componentes")
ax.legend();
En base a los valores de BIC, la configuración óptima del modelo es de entorno a 6 componentes y matriz de covarianza tipo full. Se entrena de nuevo el modelo utilizando estos valores.
# Entrenamiento del modelo GMM
# ==============================================================================
modelo_gmm = GaussianMixture(
n_components = 6,
covariance_type = 'full',
random_state = 123,
)
modelo_gmm.fit(X=datos_X)
Predicción¶
Los modelos GaussianMixture tienen varios métodos de predicción. Para utilizarlos como detectores de anomalía, interesa conocer la densidad de probabilidad que tiene cada observación acorde al modelo. Este valor puede obtenerse con el método score_samples(). En la documentación de indica que este método devuelve el logaritmo de la probabilidad, pero, estrictamente hablando, se trata del logaritmo de la densidad de probabilidad. Es por esta razón por la que, si se revierte el logaritmo con la función exponente, se obtienen valores mayores de 1.
# Predicción log probabilidad
# ==============================================================================
log_probabilidad_predicha = modelo_gmm.score_samples(X=datos_X)
log_probabilidad_predicha
# Distribución de las probabilidades (log density)
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 3.5))
sns.distplot(
log_probabilidad_predicha,
hist = False,
rug = True,
color = 'blue',
kde_kws = {'shade': True, 'linewidth': 1},
ax = ax
)
ax.set_title('Distribución predicciones')
ax.set_xlabel('Logaritmo densidad de probabilidad');
Detección de anomalías¶
Una vez calculada la densidad de probabilidad de cada observación acorde al modelo, se pueden emplear como criterio para identificar anomalías. Cabría esperar que, observaciones con una probabilidad predicha muy baja, sean anómalas.
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 probabilidades menores.
# Distribución de los valores de anomalía
# ==============================================================================
df_resultados = pd.DataFrame({
'log_probabilidad' : log_probabilidad_predicha,
'anomalia' : datos_y
})
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 3.5))
sns.boxplot(
x = 'anomalia',
y = 'log_probabilidad',
data = df_resultados,
#color = "white",
palette = 'tab10',
ax = ax
)
ax.set_title('Distribución predicciones GMM')
ax.set_ylabel('Logaritmo densidad de probabilidad')
ax.set_xlabel('clasificación (0 = normal, 1 = anomalía)');
La distribución del logaritmo de probabilidad (densidad) en el grupo de las anomalías es inferior. Sin embargo, al existir un solapamiento notable, si se clasifican las n observaciones con menor probabilidad 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 menor probabilidad predicha.
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = df_resultados \
.sort_values('log_probabilidad', ascending=True) \
.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 32% lo son realmente. El porcentaje de falsos positivos es muy elevado, el método de GMM no consigue resultados notables en este set de datos.
Una de las principales limitaciones del uso de modelos GMM como detectores de anomalías es que consideran que los datos siguen distribuciones normales multivariante. En la práctica, esta asunción es difícil que se cumpla, sobretodo a medida que aumenta el número de variables. Una estrategia que en ocasiones consigue mitigar parte de este problema es reducir la dimensionalidad de los datos con PCA y luego aplicar un modelo GMM con solo unas pocas componentes.
Reentrenamiento iterativo¶
El modelo anterior se ha entrenado empleando todas las observaciones, incluyendo potenciales anomalías. Dado que el objetivo es aprender las distribuciones de las componentes únicamente con los datos “normales”, se puede mejorar el resultado reentrenando el modelo pero excluyendo las n observaciones con menor probabilidad (potenciales anomalías).
Se repite la detección de anomalías pero, esta vez, descartando las observaciones con una densidad de probabilidad inferior al cuantil 0.01.
# Se eliminan las observaciones con una densidad de probabilidad inferior
# al cuantil 0.01
# ==============================================================================
cuantil = np.quantile(a=log_probabilidad_predicha, q=0.01)
print('Cuantil: ', cuantil)
datos_X_trimmed = datos_X.loc[log_probabilidad_predicha > cuantil, :].copy()
datos_y_trimmed = datos_y[log_probabilidad_predicha > cuantil].copy()
# Entrenamiento del modelo GMM
# ==============================================================================
modelo_gmm = GaussianMixture(
n_components = 6,
covariance_type = 'full',
random_state = 123,
)
modelo_gmm.fit(X=datos_X_trimmed)
# Matriz de confusión de la clasificación final
# ==============================================================================
df_resultados = pd.DataFrame({
'log_probabilidad': modelo_gmm.score_samples(X=datos_X),
'anomalia' : datos_y
})
df_resultados = df_resultados \
.sort_values('log_probabilidad', ascending=True) \
.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 menor probabilidad y reentrenando el modelo, 88 de las 176 observaciones identificadas como anomalías lo son realmente. Se ha reducido ligeramente el porcentaje de falsos positivos al 50%.
Ejemplo 2¶
El set de datos geyser del paquete de R MASS, contiene información sobre las erupciones del géiser Old Faithful localizado en el parque nacional de Yellowstone, Wyoming. En concreto, recoge información sobre la duración de 299 erupciones, así como el tiempo transcurrido desde la anterior.
Datos¶
# Lectura de datos
# ==============================================================================
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/' \
+ 'Estadistica-machine-learning-python/master/data/geyser.csv'
datos_X = pd.read_csv(url)
datos_X.info()
Al tratarse de datos con solo dos variables, se puede representar su distribución.
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 4))
ax.scatter(datos_X.waiting, datos_X.duration)
ax.set_title('Distribución erupciones geiser')
ax.set_ylabel('Espera')
ax.set_xlabel('Duración');
Modelo¶
En problemas con dos dimensiones (variables) como este, la identificación del número de componentes puede hacerse visualmente. Viendo los datos, parece razonable pensar que hay 3 grupos en los datos.
# Tunning del modelo GMM en base a la métrica BIC
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3.84))
n_components = range(1, 5)
covariance_types = ['spherical', 'tied', 'diag', 'full']
for covariance_type in covariance_types:
valores_bic = []
for i in n_components:
modelo = GaussianMixture(n_components=i, covariance_type=covariance_type)
modelo = modelo.fit(datos_X)
valores_bic.append(modelo.bic(datos_X))
ax.plot(n_components, valores_bic, label=covariance_type)
ax.set_title("Valores BIC")
ax.set_xlabel("Número componentes")
ax.legend();
Viendo la métrica BIC, su valor se estabiliza a partir de las 3 componentes, lo que coincide con la idea obtenida en la inspección visual de los datos.
# Entrenamiento del modelo GMM
# ==============================================================================
modelo_gmm = GaussianMixture(
n_components = 3,
covariance_type = 'diag',
random_state = 123,
)
modelo_gmm.fit(X=datos_X)
Predicción probabilidades¶
# Predicción densidad de probabilidad
# ==============================================================================
log_probabilidad_predicha = modelo_gmm.score_samples(X=datos_X)
# Mapa de densidad de probabilidad
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 4))
# Grid de valores dentro del rango observado
x = np.linspace(min(datos_X.waiting)*0.8, max(datos_X.waiting)*1.2, 1000)
y = np.linspace(min(datos_X.duration)*0.8, max(datos_X.duration)*1.2, 1000)
xx, yy = np.meshgrid(x, y)
# Densidad de probabilidad de cada valor del grid
scores = modelo_gmm.score_samples(np.c_[xx.ravel(), yy.ravel()])
scores = np.exp(scores) # Los valores están en log
ax.scatter(datos_X.waiting, datos_X.duration, alpha=0.6)
ax.contour(
xx, yy, scores.reshape(xx.shape), alpha=0.6,
levels=np.percentile(scores, np.linspace(0, 100, 10))[1:-1]
)
ax.set_title('Densidad de probabilidad del modelo');
df_resultados = datos_X.copy()
df_resultados['log_proba'] = log_probabilidad_predicha
df_resultados = df_resultados.sort_values(by='log_proba')
top_5_anomalias = df_resultados.head(10)
top_5_anomalias
# Mapa de densidad de probabilidad
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5.5, 4))
# Grid de valores dentro del rango observado
x = np.linspace(min(datos_X.waiting)*0.8, max(datos_X.waiting)*1.2, 1000)
y = np.linspace(min(datos_X.duration)*0.8, max(datos_X.duration)*1.2, 1000)
xx, yy = np.meshgrid(x, y)
# Densidad de probabilidad de cada valor del grid
scores = modelo_gmm.score_samples(np.c_[xx.ravel(), yy.ravel()])
scores = np.exp(scores) # Los valores están en log
ax.scatter(datos_X.waiting, datos_X.duration, alpha=0.6)
ax.scatter(top_5_anomalias.waiting, top_5_anomalias.duration, c="red", label='anomalía')
ax.contour(
xx, yy, scores.reshape(xx.shape),
alpha=0.6, cmap='viridis',
levels=np.percentile(scores, np.linspace(0, 100, 10))[1:-1]
)
ax.set_title('Densidad de probabilidad del modelo');
ax.legend();
Información de sesión¶
from sinfo import sinfo
sinfo()
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 Gaussian Minture Modedel (GMM) y python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py23-deteccion-anomalias-gmm-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.