Detección de Anomalías

Algoritmos de detección de anomalías en machine learning: métodos supervisados y no supervisados para identificar valores atípicos.

Palabras clave

machine learning, aprendizaje automatico, Python, algebra lineal, optimizacion, regresion lineal, clasificacion, estadistica

1 Detección de Anomalías

El objetivo de este algoritmo es evaluar un conjunto de datos normales y aprender sus patrones. Cuando llega un dato nuevo, el modelo calcula qué tan probable es que pertenezca a ese grupo común. Si esa probabilidad es extremadamente baja, el dato se etiqueta como una anomalía porque no sigue el mismo patrón o los mismos patrones que el resto de datos.

Es un algoritmo muy útil cuando muchos datos pertenecen a la clase normal y muy pocos son anomalos, por ejemplo, cuando se buscan fraudes en transacciones bancarias, la mayoría de los registros disponibles corresponderán a transacciones legítimas. El modelo no busca separar dos grupos; busca memorizar la “forma” del grupo normal y alertar si algo cae fuera.

Si utilizaramos una regresión logística por ejemplo, al tener tantos datos de una clase y pocos de la otra se sesgaria el modelo de forma natural hacia la clase habitual.

Dependiendo de cómo estén tus datos, se suelen usar tres algoritmos clásicos:

• Modelos de Densidad (Distribución Gaussiana): Estiman la probabilidad de los datos basándose en la campana de Gauss. Es el método estándar cuando las variables son continuas.
• Bosques de Aislamiento (Isolation Forests): Algoritmo basado en árboles de decisión. En lugar de aprender lo que es normal, intenta “aislar” activamente cada punto. Las anomalías requieren muy pocos pasos para ser aisladas del resto.
• One-Class SVM (Máquinas de Vector de Soporte de una Clase): Traza un límite cerrado alrededor de la zona de alta densidad de los datos normales. Cualquier punto fuera de ese límite es una anomalía.

2 Modelo de Densidad: Distribución Gaussiana

2.1 Concepto

En un modelo de densidad, el algoritmo calcula una función que dibuja dónde se concentran la mayoría de los datos (la zona de alta probabilidad). Si un dato nuevo cae en una zona donde la densidad es cercana a cero, significa que está “lejos” de lo que el modelo considera normal o probable, y por lo tanto, se dispara la alerta de anomalía.

Figura 1: Distribución de los servidores web. El algoritmo memoriza el grupo compacto (normal) y activa las alertas para los puntos que quedan fuera de la zona de alta densidad (anomalías).

Aunque intuitivamente decimos que una anomalía está “lejos”, el modelo de densidad no mide distancias geométricas como lo haría, por ejemplo, el algoritmo de los \(k\) vecinos más cercanos o k-NN. Lo que calcula es un valor de densidad de probabilidad. El modelo memoriza una distribución normal y al llegar el dato nuevo, lo introduce en la fórmula matemática y esta devuelve un número: la probabilidad \(p(X)\).

2.2 Cálculo

Para cada una de las \(n\) características, calculamos su media (\(\mu\)) y su varianza (\(\sigma^2\)) utilizando los datos de entrenamiento que normalmente no presentan anomalías:

\[\mu_j = \frac{1}{m} \sum_{i=1}^{m} x_j^{(i)}\] \[\sigma_j^2 = \frac{1}{m} \sum_{i=1}^{m} (x_j^{(i)} - \mu_j)^2\]

Para un nuevo dato \(X\) (que es un vector con \(n\) características), se calcula la probabilidad total multiplicando la probabilidad de cada una de sus partes. La ecuación de la densidad Gaussiana para una sola característica es:

\[p(x_j; \mu_j, \sigma_j^2) = \frac{1}{\sqrt{2\pi\sigma_j^2}} e^{-\frac{(x_j - \mu_j)^2}{2\sigma_j^2}}\]

Para el vector completo \(X\), la probabilidad conjunta se calcula con la productoria (\(\prod\)):

\[p(X) = \prod_{j=1}^{n} p(x_j; \mu_j, \sigma_j^2) = p(x_1) \cdot p(x_2) \cdots p(x_n)\]

Una vez que calculamos \(p(X)\), se define un umbral de tolerancia muy pequeño llamado épsilon (\(\epsilon\)):

• Si \(p(X) < \epsilon \longrightarrow\) Anomalía.
• Si \(p(X) \ge \epsilon \longrightarrow\) Dato Normal.

2.3 Ejemplo Práctico Numérico

Imagina que monitorizamos servidores web usando \(n=2\) características:

• \(x_1\): Latencia de respuesta (en ms).
• \(x_2\): Rendimiento de la CPU (en %).

Tras entrenar con servidores normales, el modelo aprende estos parámetros:

• Para la latencia (\(x_1\)): \(\mu_1 = 50\), \(\sigma_1^2 = 25\)
• Para la CPU (\(x_2\)): \(\mu_2 = 20\), \(\sigma_2^2 = 4\)

Llega un servidor con un reporte de \(X = \begin{bmatrix} 52 \\ 19 \end{bmatrix}\). Está muy cerca de las medias (50 y 20). Al calcular \(p(x_1)\) y \(p(x_2)\) con la fórmula Gaussiana, se obtiene valores relativamente altos.

\[p(X) = 0.07 \cdot 0.18 = \mathbf{0.0126}\]

Si el umbral es \(\epsilon = 0.001\), como \(0.0126 > 0.001\), el sistema determina que el servidor está sano.

Llega un nuevo servidor con \(X = \begin{bmatrix} 95 \\ 85 \end{bmatrix}\) (alta latencia y CPU al límite).Como 95 está lejos de la media de latencia (50) y 85 está lejos de la CPU promedio (20), el exponente negativo de la fórmula de Gauss se dispara hacia abajo, dando probabilidades muy pequeñas.

\[p(X) = (1.2 \times 10^{-7}) \cdot (3.1 \times 10^{-15}) = \mathbf{3.72 \times 10^{-22}}\]

Como este valor es críticamente menor que \(\epsilon = 0.001\), el algoritmo lo clasifica como anomalía.

2.4 Desvanecimiento del producto

Si cada probabilidad individual \(p(x_j)\) es un número entre 0 y 1 (por ejemplo, \(0.1\)), al multiplicar 2 características tienes \(0.1 \times 0.1 = 0.01\). Pero si tienes 50 características, \(0.1^{50}\) se convierte en un número muy cercano a cero (\(10^{-50}\)), no porque el dato sea una anomalía, sino por multiplicación matemática.

Para evitar que las probabilidades colapsen hacia cero, en la práctica nunca se calcula el producto directo. Se aplica el logaritmo natural a la función de densidad.Por propiedades de los logaritmos, una multiplicación se transforma en una suma:

\[\log p(X) = \log(p(x_1) \cdot p(x_2) \cdots p(x_n)) = \sum_{j=1}^{n} \log p(x_j)\]

Al transformar el producto en una suma, los números ya no se desvanecen exponencialmente. En lugar de buscar un \(\epsilon\) ridículamente pequeño (como \(10^{-50}\)), buscas un umbral de coste logarítmico negativo (por ejemplo, \(\epsilon = -25\)).

Imagina que estamos analizando un sistema con 5 características (\(n=5\)). Tras calcular la probabilidad de cada una de ellas para un dato sospechoso, la fórmula de Gauss nos devuelve las siguientes densidades de probabilidad individuales:

• \(p(x_1) = 0.01\)
• \(p(x_2) = 0.03\)
• \(p(x_3) = 0.005\)
• \(p(x_4) = 0.02\)
• \(p(x_5) = 0.01\)

Calculamos los logaritmos de nuestras 5 probabilidades (redondeando a dos decimales):

• \(\ln(0.01) = -4.61\)
• \(\ln(0.03) = -3.51\)
• \(\ln(0.005) = -5.30\)
• \(\ln(0.02) = -3.91\)
• \(\ln(0.01) = -4.61\)

Ahora, la propiedad de los logaritmos nos permite sumar estos valores en lugar de multiplicarlos:

\[\ln p(X) = \sum_{j=1}^{5} \ln p(x_j)\] \[\ln p(X) = (-4.61) + (-3.51) + (-5.30) + (-3.91) + (-4.61) = \mathbf{-21.94}\]

Ahora debemos comprar el resultado con un épsilon en formato de puntuación negativa, por ejemplo: \(\ln \epsilon = -25\).

Como el coste logarítmico de nuestro dato (\(-21.94\)) es mayor (menos negativo) que el umbral de peligro (\(-25\)), el algoritmo determina que el dato sigue estando dentro del rango normal tolerado.

2.5 Entrenamiento

2.5.1 Paso 1. Obtener los datos de la distribución de ejemplos sanos

Recogemos datos de 5 turbinas que sabemos que funcionan a la perfección (\(m=5\)). Como es entrenamiento para el modelo de densidad, aquí no hay etiquetas.La matriz de datos \(X\) (orientada por columnas, donde cada columna es una turbina y cada fila una característica) es:

\[X = \begin{bmatrix} 48 & 50 & 52 & 49 & 51 \\ 12 & 10 & 11 & 13 & 9 \end{bmatrix}\]

Calcular el promedio para cada fila (característica):

\[\mu_1 = \frac{48 + 50 + 52 + 49 + 51}{5} = \frac{250}{5} = \mathbf{50}\]

\[\mu_2 = \frac{12 + 10 + 11 + 13 + 9}{5} = \frac{55}{5} = \mathbf{11}\]

Calcular la varianza:

\[\sigma_1^2 = \frac{(48-50)^2 + (50-50)^2 + (52-50)^2 + (49-50)^2 + (51-50)^2}{5} = \frac{4 + 0 + 4 + 1 + 1}{5} = \frac{10}{5} = \mathbf{2}\] \[\sigma_2^2 = \frac{(12-11)^2 + (10-11)^2 + (11-11)^2 + (13-11)^2 + (9-11)^2}{5} = \frac{1 + 1 + 0 + 4 + 4}{5} = \frac{10}{5} = \mathbf{2}\]

• Temperatura: \(\mu_1 = 50, \sigma_1^2 = 2\)
• Vibración: \(\mu_2 = 11, \sigma_2^2 = 2\)

2.5.2 Paso 2. El Set de Validación

Para encontrar el umbral de peligro, hay que probar el modelo en un conjunto de validación que sí esté etiquetado. Supongamos que tenemos 3 casos del pasado:

• Caso A (Turbina normal): \(X_A = \begin{bmatrix} 50 \\ 12 \end{bmatrix}\) con etiqueta \(Y_A = 0\) (Normal)
• Caso B (Anomalía - Sobrecalentamiento): \(X_B = \begin{bmatrix} 56 \\ 11 \end{bmatrix}\) con etiqueta \(Y_B = 1\) (Anomalía)
• Caso C (Anomalía - Vibración extrema): \(X_C = \begin{bmatrix} 49 \\ 17 \end{bmatrix}\) con etiqueta \(Y_C = 1\) (Anomalía)

La fórmula del logaritmo de una densidad Gaussiana es:

\[\ln p(x_j) = \ln\left(\frac{1}{\sqrt{2\pi\sigma_j^2}}\right) - \frac{(x_j - \mu_j)^2}{2\sigma_j^2}\]

Como en este caso ambas varianzas valen \(\sigma^2 = 2\), el primer bloque de la ecuación siempre valdrá:

\[\ln\left(\frac{1}{\sqrt{2\pi \cdot 2}}\right) = \ln(0.282) = \mathbf{-1.265}\]

Por lo tanto, la fórmula simplificada para este problema exacto es:

\[\ln p(x_j) = -1.265 - \frac{(x_j - \mu_j)^2}{4}\]

2.5.3 Paso 3: Calcular las puntuaciones de Validación

Ahora disponemos de tres casos, uno que es normal y dos que presentan anomalía. Al ser dos caracteristicas aplicamos la fórmula a cada uno de los casos \(\ln p(X) = \ln p(x_1) + \ln p(x_2)\).

Evaluando el Caso A:

• \(X_A = \begin{bmatrix} 50 \\ 12 \end{bmatrix}\)
• \(\ln p(x_1) = -1.265 - \frac{(50-50)^2}{4} = -1.265 - 0 = -1.265\)
• \(\ln p(x_2) = -1.265 - \frac{(12-11)^2}{4} = -1.265 - 0.25 = -1.515\)
• Puntuación Total Caso A: \(-1.265 + (-1.515) = \mathbf{-2.78}\)

Evaluando el Caso B:

• \(X_B = \begin{bmatrix} 56 \\ 11 \end{bmatrix}\)
• \(\ln p(x_1) = -1.265 - \frac{(56-50)^2}{4} = -1.265 - \frac{36}{4} = -1.265 - 9 = -10.265\)
• \(\ln p(x_2) = -1.265 - \frac{(11-11)^2}{4} = -1.265 - 0 = -1.265\)
• Puntuación Total Caso B: \(-10.265 + (-1.265) = \mathbf{-11.53}\)

Evaluando el Caso C:

• \(X_C = \begin{bmatrix} 49 \\ 17 \end{bmatrix}\)
• \(\ln p(x_1) = -1.265 - \frac{(49-50)^2}{4} = -1.265 - 0.25 = -1.515\)
• \(\ln p(x_2) = -1.265 - \frac{(17-11)^2}{4} = -1.265 - \frac{36}{4} = -10.265\)
• Puntuación Total Caso C: \(-1.515 + (-10.265) = \mathbf{-11.78}\)

Como vemos a más bajo sea el valor menor es la probabilidad de que sea un dato normal. Y será más bajo cuanto menos comunes sean las propiedades del vector.

2.5.4 Paso 4: Encontrar el Umbral

Ahora en base a los dato de validación y sus puntuaciones necesitamos encontrar un umbral que sea lo más optimos posible, para ello se utiliza el F1 Score a partir de las puntuaciones de validación:

• Caso A (Real: 0): Nota = \(-2.78\)
• Caso B (Real: 1): Nota = \(-11.53\)
• Caso C (Real: 1): Nota = \(-11.78\)

Se elije un humbral y por cada umbral se cálcula:

\[F_1 = 2 \cdot \frac{\text{Precisión} \cdot \text{Exhaustividad}}{\text{Precisión} + \text{Exhaustividad}}\]

\[\text{Precisión} = \frac{\text{VP}}{\text{VP} + \text{FP}} \] \[\text{Exhaustividad (Recall)} = \frac{\text{VP}}{\text{VP} + \text{FN}}\]

• VP (Verdaderos Positivos): Casos en los que el modelo predijo que había una anomalía/fraude (\(1\)) y la realidad demostró que sí era una anomalía/fraude (\(1\)). Es decir, alertas que resultaron ser totalmente ciertas.
• FP (Falsos Positivos): Casos en los que el modelo predijo que había una anomalía/fraude (\(1\)) pero en la realidad el dato estaba completamente sano/normal (\(0\)). Esto representa las falsas alarmas.
• FN (Falso Negativo): Desastre silencioso. El modelo dice “Sano” pero hay peligro. (Te cuesta dinero o seguridad).

El F1-Score devolverá un valor entre 0 y 1

• Si da 1: Significa que tu umbral logró un modelo perfecto (cero falsas alarmas y cero fallos escapados).
• Si da 0: El modelo es completamente inútil.

Vamos a seleccionar tres umbrales y a probar uno a uno. Si la nota es menor (más negativa) que el umbral, el modelo predice anomalía (\(1\)). Si es mayor o igual, predice normal (\(0\)). A más pequeño menor es la probabilidad de pertenecer a los datos normales.

2.5.4.1 \(\ln \epsilon = -12\)

Predicciones:

• Caso A (\(-2.78 \ge -12\)) \(\rightarrow\) Predice: \(0\) (Normal) Correcto
• Caso B (\(-11.53 \ge -12\)) \(\rightarrow\) Predice: \(0\) (Normal) Error
• Caso C (\(-11.78 \ge -12\)) \(\rightarrow\) Predice: \(0\) (Normal) Error

Recuento:

• \(\text{VP} = 0\)
• \(\text{FP} = 0\)
• \(\text{FN} = 2\)

Metricas:

• Precisión: Indefinida (\(0/0\)), porque no envió ninguna alarma.
• Exhaustividad (Recall): \(\frac{0}{0 + 2} = \mathbf{0.0}\) (No atrapó nada).
• F1-Score: \(\mathbf{0.0}\)

2.5.4.2 \(\ln \epsilon = -11.6\)

Predicciones:

• Caso A (\(-2.78 \ge -11.6\)) \(\rightarrow\) Predice: \(0\) (Normal) Correcto
• Caso B (\(-11.53 \ge -11.6\)) \(\rightarrow\) Predice: \(0\) (Normal) Error
• Caso C (\(-11.78 < -11.6\)) \(\rightarrow\) Predice: \(1\) (Anomalía) Correcto

Recuento:

• \(\text{VP} = 1\)
• \(\text{FP} = 0\)
• \(\text{FN} = 1\) El caso B.

Metricas:

• Precisión: $ = $1
• Exhaustividad (Recall): \(\frac{1}{1 + 1} = \mathbf{0.5}\)
• F1-Score: \(2 \cdot \frac{1.0 \cdot 0.5}{1.0 + 0.5} = \mathbf{0.66}\)

2.5.4.3 \(\ln \epsilon = -10\)

Predicciones:

• Caso A (\(-2.78 \ge -10\)) \(\rightarrow\) Predice: \(0\) (Normal) Correcto
• Caso B (\(-11.53 < -10\)) \(\rightarrow\) Predice: \(1\) (Anomalía) Correcto
• Caso C (\(-11.78 < -10\)) \(\rightarrow\) Predice: \(1\) (Anomalía) Correcto

Recuento:

• \(\text{VP} = 2\)
• \(\text{FP} = 0\)
• \(\text{FN} = 0\)

Metricas:

• Precisión: \(\frac{2}{2 + 0} = \mathbf{1.0}\)
• Exhaustividad (Recall): \(\frac{2}{2 + 0} = \mathbf{1.0}\)
• F1-Score: \(2 \cdot \frac{1.0 \cdot 1.0}{1.0 + 1.0} = \mathbf{1.0}\)

El script automático se quedaría con el umbral de \(-10\) porque su F1-Score es de 1.0. Gráficamente, el \(-10\) actúa como una línea divisoria perfecta que deja al caso sano a un lado y a los dos casos rotos al otro, protegiendo la planta nuclear sin molestar a los ingenieros con falsas alertas.