Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Más sobre forecasting en: cienciadedatos.net
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.Las librerías utilizadas en este documento son:
# 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__)
# Descarga de datos
# ==============================================================================
data = fetch_dataset("ett_m1")
print(f"Data shape: {data.shape}")
data.head(3)
# 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.
# 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)
# 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();
# 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)}")
# 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();
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.
# 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))
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.
# 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
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.
# 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)
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.
# 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()
display(top_residuals_train)
top_residuals_train["is_anomaly"].value_counts()
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}%")
El mismo proceso se repite para el conjunto de test, donde hay un total de 4 anomalías.
# 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)
# 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()
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()
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.import session_info
session_info.show(html=False)
¿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} }
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.