Detección de anomalías en series temporales

Si te gusta  Skforecast , ayúdanos dándonos una estrella en  GitHub! ⭐️

Detección de anomalías en series temporales

Joaquin Amat Rodrigo, Javier Escobar Ortiz
Diciembre, 2024

Introducción

La detección de anomalías implica identificar valores, patrones o eventos que se desvían del comportamiento esperado. En el contexto de las series temporales, estos valores atípicos pueden tener diversas causas, como cambios en el proceso subyacente de generación de datos, errores de medición o acontecimientos externos inesperados. Detectar anomalías en las series temporales es especialmente difícil debido a las características inherentes a estos datos, como las dependencias temporales, las tendencias, la estacionalidad y el ruido. Estos factores pueden dificultar la distinción entre anomalías y variaciones normales.

A pesar de estas dificultades, la detección de anomalías en las series temporales es fundamental. Si no se detectan, las anomalías pueden afectar significativamente al rendimiento de los modelos de forecasting, sesgando sus predicciones y comprometiendo su precisión. Por lo tanto, es esencial identificar y analizar cuidadosamente las anomalías para mantener la fiabilidad de los modelos y mejorar la calidad de las percepciones derivadas de los datos.

Este artículo explora un enfoque específico para la detección de anomalías en las series temporales: la medición de los residuos de un modelo de forecasting. Los residuos representan las diferencias entre los valores observados y las correspondientes predicciones del modelo, y ofrecen una perspectiva valiosa para identificar valores anómalos que el modelo tiene dificultades para predecir con exactitud.

✎ Note

Existen numerosos enfoques para la detección de anomalías en las series temporales, incluidos los métodos estadísticos, la descomposición de tendencias estacionales, las técnicas de aprendizaje automático como la clasificación y la agrupación, así como los autocodificadores, entre otros. Cada método tiene sus puntos fuertes y sus limitaciones, y la elección del enfoque debe ajustarse a las características específicas de los datos y a los objetivos del análisis.

Librerías

Las librerías utilizadas en este documento son:

In [20]:
# Procesado de datos
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme

# Forecastig
# ==============================================================================
import skforecast
import sklearn
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster
from skforecast.preprocessing import RollingFeatures

# Warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

print('Versión skforecast:', skforecast.__version__)
print('Versión sklearn:', sklearn.__version__)
Versión skforecast: 0.15.0
Versión sklearn: 1.5.1

Datos

In [21]:
# Descarga de datos
# ==============================================================================
data = fetch_dataset("ett_m1")
print(f"Data shape: {data.shape}")
data.head(3)
ett_m1
------
Data from an electricity transformer station was collected between July 2016 and
July 2018 (2 years × 365 days × 24 hours × 4 intervals per hour = 70,080 data
points). Each data point consists of 8 features, including the date of the
point, the predictive value "Oil Temperature (OT)", and 6 different types of
external power load features: High UseFul Load (HUFL), High UseLess Load (HULL),
Middle UseFul Load (MUFL), Middle UseLess Load (MULL), Low UseFul Load (LUFL),
Low UseLess Load(LULL).
Zhou, Haoyi & Zhang, Shanghang & Peng, Jieqi & Zhang, Shuai & Li, Jianxin &
Xiong, Hui & Zhang, Wancai. (2020). Informer: Beyond Efficient Transformer for
Long Sequence Time-Series Forecasting.
[10.48550/arXiv.2012.07436](https://arxiv.org/abs/2012.07436).
https://github.com/zhouhaoyi/ETDataset
Shape of the dataset: (69680, 7)
Data shape: (69680, 7)
Out[21]:
HUFL HULL MUFL MULL LUFL LULL OT
date
2016-07-01 00:00:00 5.827 2.009 1.599 0.462 4.203 1.340 30.531000
2016-07-01 00:15:00 5.760 2.076 1.492 0.426 4.264 1.401 30.459999
2016-07-01 00:30:00 5.760 1.942 1.492 0.391 4.234 1.310 30.038000
In [22]:
# Selección de datos
# ==============================================================================
data = data.loc["20170101":, ["OT"]]

Se añaden algunas anomalías sintéticas a los datos de la serie temporal para ilustrar el proceso de detección.

In [23]:
# Anomalías sintéticas
# ==============================================================================
rng = np.random.default_rng(4631789)
n_anomalies = 20
anomalies_loc = data.index[rng.uniform(low=0, high=len(data), size=n_anomalies).astype(int)]
data.loc[anomalies_loc, "OT"] = data.loc[anomalies_loc, "OT"] * rng.choice([0.5, 2], size=n_anomalies)
In [24]:
# Gráfico
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(data["OT"])
ax.scatter(anomalies_loc, data.loc[anomalies_loc, "OT"], color='orange', label='Anomalies', zorder=2)
ax.set_title('OT')
ax.legend();
In [25]:
# Partición de datos
# ==============================================================================
end_train = '2018-04-03 23:59:00'
data_train = data.loc[: end_train, :]
data_test  = data.loc[end_train:, :]

anomalies_loc_train = [date for date in anomalies_loc if date <= pd.to_datetime(end_train)]
anomalies_loc_test  = [date for date in anomalies_loc if date > pd.to_datetime(end_train)]

print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"    Anomalies train: {len(anomalies_loc_train)}")
print(f"Dates test  : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")
print(f"    Anomalies test : {len(anomalies_loc_test)}")
Dates train : 2017-01-01 00:00:00 --- 2018-04-03 23:45:00  (n=43968)
    Anomalies train: 16
Dates test  : 2018-04-04 00:00:00 --- 2018-06-26 19:45:00  (n=8048)
    Anomalies test : 4
In [26]:
# Plot partitions
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(data_train['OT'], label='Train')
ax.plot(data_test['OT'], label='Test')
ax.scatter(anomalies_loc, data.loc[anomalies_loc, "OT"], color='orange', label='Anomalies', zorder=3)
ax.set_title('OT')
ax.legend();

Feature engineering

Se crean características exógenas adicionales basadas en:

  • Valores pasados de la serie exógena.

  • Estadísticos móviles de la serie exógena.

  • Variables de calendario.

In [27]:
# Variables de calendario
# ==============================================================================
features_to_extract = [
    'year',
    'month',
    'week',
    'day_of_week',
    'hour'
]
calendar_transformer = DatetimeFeatures(
    variables           = 'index',
    features_to_extract = features_to_extract,
    drop_original       = False,
)

# Codificación cíclica
# ==============================================================================
features_to_encode = [
    "month",
    "week",
    "day_of_week",
    "hour",
]
max_values = {
    "month": 12,
    "week": 52,
    "day_of_week": 7,
    "hour": 24,
}
cyclical_encoder = CyclicalFeatures(
                        variables     = features_to_encode,
                        max_values    = max_values,
                        drop_original = True
                   )

exog_transformer = make_pipeline(
                        calendar_transformer,
                        cyclical_encoder
                   )
display(exog_transformer)

data = exog_transformer.fit_transform(data)
# Remove rows with NaNs created by lag features
data = data.dropna()
exog_features = data.columns.difference(['OT']).tolist()
display(data.head(3))
Pipeline(steps=[('datetimefeatures',
                 DatetimeFeatures(drop_original=False,
                                  features_to_extract=['year', 'month', 'week',
                                                       'day_of_week', 'hour'],
                                  variables='index')),
                ('cyclicalfeatures',
                 CyclicalFeatures(drop_original=True,
                                  max_values={'day_of_week': 7, 'hour': 24,
                                              'month': 12, 'week': 52},
                                  variables=['month', 'week', 'day_of_week',
                                             'hour']))])
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
OT year month_sin month_cos week_sin week_cos day_of_week_sin day_of_week_cos hour_sin hour_cos
date
2017-01-01 00:00:00 10.200 2017 0.5 0.866025 -2.449294e-16 1.0 -0.781831 0.62349 0.0 1.0
2017-01-01 00:15:00 9.989 2017 0.5 0.866025 -2.449294e-16 1.0 -0.781831 0.62349 0.0 1.0
2017-01-01 00:30:00 10.130 2017 0.5 0.866025 -2.449294e-16 1.0 -0.781831 0.62349 0.0 1.0

Entrenamiento del modelo

Se entrena un forecaster para que sus predicciones puedan compararse con los valores observados. Cuanto más eficazmente aprenda el modelo los patrones subyacentes en los datos, más precisas serán sus predicciones. Esta mayor precisión aumenta la confianza en que las desviaciones entre los valores observados y los previstos representan auténticas anomalías.

In [28]:
# Forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=24)
forecaster = ForecasterRecursive(
                 regressor       = HistGradientBoostingRegressor(random_state=15926, loss='quantile', quantile=0.5),
                 lags            = [1, 3, 4, 5, 6, 13, 14, 15, 16, 17],
                 differentiation = 1,
                 transformer_y   = StandardScaler(),
                 window_features = window_features,
             )
forecaster.fit(
    y    = data.loc[:end_train, 'OT'],
    exog = data.loc[:end_train, exog_features]
)
forecaster
Out[28]:

ForecasterRecursive

General Information
  • Regressor: HistGradientBoostingRegressor
  • Lags: [ 1 3 4 5 6 13 14 15 16 17]
  • Window features: ['roll_mean_24', 'roll_min_24', 'roll_max_24']
  • Window size: 25
  • Exogenous included: True
  • Weight function included: False
  • Differentiation order: 1
  • Creation date: 2025-02-03 16:41:14
  • Last fit date: 2025-02-03 16:41:15
  • Skforecast version: 0.15.0
  • Python version: 3.12.5
  • Forecaster id: None
Exogenous Variables
    day_of_week_cos, day_of_week_sin, hour_cos, hour_sin, month_cos, month_sin, week_cos, week_sin, year
Data Transformations
  • Transformer for y: StandardScaler()
  • Transformer for exog: None
Training Information
  • Training range: [Timestamp('2017-01-01 00:00:00'), Timestamp('2018-04-03 23:45:00')]
  • Training index type: DatetimeIndex
  • Training index frequency: 15min
Regressor Parameters
    {'categorical_features': 'warn', 'early_stopping': 'auto', 'interaction_cst': None, 'l2_regularization': 0.0, 'learning_rate': 0.1, 'loss': 'quantile', 'max_bins': 255, 'max_depth': None, 'max_features': 1.0, 'max_iter': 100, 'max_leaf_nodes': 31, 'min_samples_leaf': 20, 'monotonic_cst': None, 'n_iter_no_change': 10, 'quantile': 0.5, 'random_state': 15926, 'scoring': 'loss', 'tol': 1e-07, 'validation_fraction': 0.1, 'verbose': 0, 'warm_start': False}
Fit Kwargs
    {}

🛈 API Reference    🗎 User Guide

Detección de anomalías

Una vez entrenado el modelo, se generan predicciones tanto para el conjunto de entrenamiento como para el conjunto de prueba. Esto se consigue mediante un proceso de backtesting. Para este ejemplo se supone que se predicen 24 pasos adelante en cada iteración.

In [30]:
# Backtesting con datos de entrenamiento
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = None, # Utilizar todos los datos asumiendo que el modelo ya está entrenado
        steps              = 24,   # Todas las horas del siguiente día
        differentiation    = 1,
     )

metric_train, predictions_train = backtesting_forecaster(
                                    forecaster    = forecaster,
                                    y             = data.loc[:end_train, 'OT'],
                                    exog          = data.loc[:end_train, exog_features],
                                    cv            = cv,
                                    metric        = 'mean_absolute_error',
                                    n_jobs        = 'auto',
                                    verbose       = False,
                                    show_progress = True
                                )

results_train = predictions_train.copy()
results_train['y_true'] = data.loc[:end_train, 'OT']
results_train['residual'] = results_train['y_true'] - results_train['pred']
results_train['residual_abs'] = np.abs(results_train['residual'])
results_train.head(3)
Out[30]:
pred y_true residual residual_abs
2017-01-01 06:15:00 10.101177 10.833 0.731823 0.731823
2017-01-01 06:30:00 10.094863 11.537 1.442136 1.442136
2017-01-01 06:45:00 10.089609 11.537 1.447390 1.447390

En el conjunto de entrenamiento, hay un total de 16 anomalías. A continuación, se identifican los 16 residuos más grandes y se extraen los puntos temporales correspondientes. Luego, estos puntos temporales se comparan con las anomalías reales para determinar la capacidad del modelo para detectarlas.

In [31]:
# Gráfico predicciones, valor real y residuos
# ==============================================================================
fig, axs = plt.subplots(
    2, 1, figsize=(7, 5), gridspec_kw={"height_ratios": [2, 1]}, sharex=True
)
axs[0].plot(data.loc[:end_train, "OT"], label="Real value")
axs[0].plot(predictions_train["pred"], label="Prediction")
axs[0].scatter(
    anomalies_loc_train,
    data.loc[anomalies_loc_train, "OT"],
    zorder=2,
    color="orange",
    label="Anomalies",
)
axs[0].set_title("Backtesting on train data")
axs[0].legend()
axs[1].plot(results_train["residual_abs"], label="Absolute residuals", zorder=2)
axs[1].set_title("Absolute residuals")
top_residuals_train = (
    results_train.sort_values("residual_abs", ascending=False).head(len(anomalies_loc_train))
)
top_residuals_train["is_anomaly"] = top_residuals_train.index.isin(anomalies_loc_train)
axs[1].scatter(
    anomalies_loc_train,
    results_train.loc[anomalies_loc_train, "residual_abs"],
    zorder=2,
    color="orange",
    label="Anomalies",
)
plt.tight_layout()
In [32]:
display(top_residuals_train)
top_residuals_train["is_anomaly"].value_counts()
pred y_true residual residual_abs is_anomaly
2017-06-21 19:15:00 19.676958 38.549999 18.873041 18.873041 True
2017-08-23 03:00:00 15.963411 32.360001 16.396589 16.396589 True
2017-04-07 21:00:00 13.038298 28.420000 15.381702 15.381702 True
2017-06-05 00:45:00 14.800458 28.982000 14.181542 14.181542 True
2017-03-04 17:15:00 13.874956 26.872000 12.997044 12.997044 True
2017-06-01 22:30:00 22.993038 10.482000 -12.511038 12.511038 False
2017-03-05 20:00:00 12.504316 24.340000 11.835685 11.835685 True
2017-06-01 22:45:00 22.969251 11.396000 -11.573251 11.573251 False
2017-05-12 04:30:00 17.299185 7.105000 -10.194185 10.194185 False
2017-07-13 22:45:00 21.299531 11.115000 -10.184532 10.184532 True
2018-03-15 21:00:00 11.082166 0.985000 -10.097166 10.097166 False
2017-11-10 02:45:00 8.333463 18.430000 10.096537 10.096537 True
2017-04-17 04:30:00 17.505625 7.457000 -10.048625 10.048625 False
2018-03-15 21:15:00 11.061701 1.126000 -9.935701 9.935701 False
2018-03-15 20:45:00 11.108253 1.196000 -9.912253 9.912253 False
2017-05-11 17:15:00 23.697282 13.788000 -9.909282 9.909282 False
Out[32]:
is_anomaly
True     8
False    8
Name: count, dtype: int64
In [33]:
print(f"Percentage of true positives: {100 * top_residuals_train['is_anomaly'].mean():.2f}%")
print(f"Percentage of false positives: {100 * (~top_residuals_train['is_anomaly']).mean():.2f}%")
Percentage of true positives: 50.00%
Percentage of false positives: 50.00%

El mismo proceso se repite para el conjunto de test, donde hay un total de 4 anomalías.

In [34]:
# Backtesting en test
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = len(data_train),
        steps              = 24,  # all hours of next day
        refit              = False,
        differentiation    = 1,
     )

metric_test, predictions_test = backtesting_forecaster(
                                    forecaster    = forecaster,
                                    y             = data.loc[:, 'OT'],
                                    exog          = data.loc[:, exog_features],
                                    cv            = cv,
                                    metric        = 'mean_absolute_error',
                                    n_jobs        = 'auto',
                                    verbose       = False,
                                    show_progress = True
                                )
results_test = predictions_test.copy()
results_test['y_true'] = data.loc[end_train:, 'OT']
results_test['residual'] = results_test['y_true'] - results_test['pred']
results_test['residual_abs'] = np.abs(results_test['residual'])
results_test.head(3)
Out[34]:
pred y_true residual residual_abs
2018-04-04 00:00:00 10.539706 9.919 -0.620706 0.620706
2018-04-04 00:15:00 10.515761 9.004 -1.511761 1.511761
2018-04-04 00:30:00 10.487109 8.723 -1.764109 1.764109
In [35]:
# Gráfico predicciones, valor real y residuos
# ==============================================================================
fig, axs = plt.subplots(
    2, 1, figsize=(8, 5), gridspec_kw={"height_ratios": [2, 1]}, sharex=True
)
axs[0].plot(data.loc[end_train:, "OT"], label="Real value")
axs[0].plot(results_test["pred"], label="Prediction")
axs[0].scatter(
    anomalies_loc_test,
    data.loc[anomalies_loc_test, "OT"],
    color="orange",
    zorder=2,
    label="Anomalies",
)
axs[0].set_title("Backtesting on test data")
axs[0].legend()
axs[1].plot(results_test["residual_abs"], label="Absolute residuals")
axs[1].set_title("Absolute residuals")
axs[1].scatter(
    anomalies_loc_test,
    results_test.loc[anomalies_loc_test, "residual_abs"],
    color="orange",
    zorder=2,
    label="Anomalies",
)
plt.tight_layout()
In [36]:
top_residuals_test = results_test.sort_values('residual_abs', ascending=False).head(len(anomalies_loc_test))
top_residuals_test['is_anomaly'] = top_residuals_test.index.isin(anomalies_loc_test)
display(top_residuals_test)
top_residuals_test['is_anomaly'].value_counts()
pred y_true residual residual_abs is_anomaly
2018-05-07 20:45:00 6.913693 17.728001 10.814307 10.814307 True
2018-05-18 16:15:00 14.043144 5.557000 -8.486144 8.486144 False
2018-05-18 16:30:00 14.080936 6.402000 -7.678936 7.678936 False
2018-05-26 02:15:00 7.083487 14.632000 7.548513 7.548513 True
Out[36]:
is_anomaly
True     2
False    2
Name: count, dtype: int64

Tanto en el conjunto de como en el de test, las observaciones con mayores residuos son anomalías verdaderas. Sin embargo, al seleccionar tantos residuos superiores como anomalías hay en los datos, sólo el 50% de ellos son anomalías verdaderas, siendo el resto falsos positivos.

Este enfoque tiene una limitación importante: las anomalías en los datos pueden distorsionar las predicciones si se incluyen en los lags utilizados por el modelo. A menudo se producen tantos falsos positivos como anomalías verdaderas, ya que la precisión del modelo se ve comprometida por las propias anomalías. Una forma de mitigar este problema es eliminar/reemplazar iterativamente las anomalías detectadas en los datos y volver a entrenar el modelo. Este proceso puede repetirse hasta que el modelo deje de detectar anomalías, momento en el que es probable que las anomalías restantes sean reales.

Warning

La eliminación de anomalías debe hacerse con precaución. Dado que se han identificado porque el modelo es incapaz de predecirlas con exactitud, es probable que su eliminación mejore el rendimiento del modelo. Sin embargo, si las anomalías no se deben a errores o problemas de medición, eliminarlas podría suponer la pérdida de información valiosa. Además, si se prevé que en el futuro se produzcan anomalías similares, el modelo no podrá aprender de ellas y su rendimiento se verá comprometido.

Información de sesión

In [37]:
import session_info
session_info.show(html=False)
-----
feature_engine      1.8.1
matplotlib          3.9.2
numpy               2.0.0
pandas              2.2.3
session_info        1.0.0
skforecast          0.15.0
sklearn             1.5.1
-----
IPython             8.27.0
jupyter_client      8.6.3
jupyter_core        5.7.2
notebook            6.4.12
-----
Python 3.12.5 | packaged by Anaconda, Inc. | (main, Sep 12 2024, 18:27:27) [GCC 11.2.0]
Linux-5.15.0-1075-aws-x86_64-with-glibc2.31
-----
Session information updated at 2025-02-03 16:42

Instrucciones para citar

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Detección de anomalías en series temporales por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible bajo una licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://cienciadedatos.net/documentos/py62-deteccion-anomalias-series-temporales.html

¿Cómo citar skforecast?

Si utilizas skforecast, te agradeceríamos mucho que lo cites. ¡Muchas gracias!

Zenodo:

Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2024). skforecast (v0.14.0). Zenodo. https://doi.org/10.5281/zenodo.8382788

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.14.0) [Computer software]. https://doi.org/10.5281/zenodo.8382788

BibTeX:

@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.14.0}, month = {11}, year = {2024}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }

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

Tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Become a GitHub Sponsor Become a GitHub Sponsor

Creative Commons Licence

Este documento creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz 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.

  • No-Comercial: No puedes utilizar el material para fines comerciales.

  • Compartir-Igual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.