Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Forecasting probabilístico con machine learning
Más sobre forecasting en: cienciadedatos.net
- Skforecast: forecasting series temporales con machine learning
- Forecasting series temporales con gradient boosting: XGBoost, LightGBM y CatBoost
- Global Forecasting: Multi-series forecasting
- Predicción de la demanda eléctrica con machine learning
- Modelos ARIMA y SARIMAX
- Demanda intermitente con skforecast
Introducción¶
Al tratar de anticipar valores futuros de una serie temporal, la mayoría de los modelos de forecasting intentan predecir cuál será el valor más probable, esto se llama point-forecasting. Aunque conocer de antemano el valor esperado de una serie temporal es útil en casi todos los casos de negocio, este tipo de predicción no proporciona información sobre la confianza del modelo ni sobre la incertidumbre de sus predicciones.
El forecasting probabilístico, a diferencia del point-forecasting, es una familia de técnicas que permiten predecir la distribución esperada de la variable respuesta en lugar de un único valor puntual. Este tipo de forecasting proporciona información muy valiosa ya que permite crear intervalos de predicción, es decir, el rango de valores donde es más probable que pueda estar el valor real. Más formalmente, un intervalo de predicción define el intervalo dentro del cual se espera encontrar el verdadero valor de la variable respuesta con una determinada probabilidad.
Existen múltiples formas de estimar intervalos de predicción en modelos de forecasting, la mayoría de las cuales requieren que los residuos (errores) del modelo sigan una distribución normal. Cuando no se puede asumir esta propiedad, dos alternativas comúnmente utilizadas son el bootstrapping y la regresión cuantílica. Para ilustrar cómo la librería skforecast permite estimar intervalos de predicción en modelos de forecasting multi-step se muestran dos ejemplo:
Intervalos de predicción basados en bootstrapping para un modelo de tipo recursive-multi-step forecaster.
Intervalos de predicción basados en regresión cuantílica para un modelo de tipo direct-multi-step forecaster.
⚠ Warning
Tal y como describe Rob J Hyndman en su blog, en los casos reales, casi todos los intervalos de predicción resultan ser demasiado estrechos. Por ejemplo, intervalos teóricos del 95% solo suelen conseguir una cobertura real de entre el 71% y el 87%. Este fenómeno surge debido a que estos intervalos no contemplan todas las fuentes de incertidumbre que, en el caso de modelos de *forecasting*, suelen ser de 4 tipos:- El término de error aleatorio
- Las estimación de parámetros
- La elección del modelo
- El proceso de predicción de valores futuros
✎ Note
Si es la primera vez que usa skforecast, visite Skforecast: forecasting de series temporales con Python y Scikit-learn para una introducción rápida.✎ Note
Conformal prediction es un método relativamente nuevo que permite estimar la incertidumbre aociada a las predicciones de modelos de machine learning. Este método está en la hoja de ruta de skforecast, pero aún no está disponible.Intervalos de predicción utilizando bootstrapping de los residuos¶
El error en la predicción del siguiente valor de una serie (one-step-ahead forecast) se define como $e_t = y_t - \hat{y}_{t|t-1}$. Asumiendo que los errores futuros serán similares a los errores pasados, es posible simular diferentes predicciones tomando muestras de los errores vistos previamente en el pasado (es decir, los residuos) y agregándolos a las predicciones.
Al hacer esto repetidamente, se crea una colección de predicciones ligeramente diferentes (posibles caminos futuros), que representan la varianza esperada en el proceso de forecasting.
Finalmente, los intervalos de predicción se crean calculando los percentiles $\alpha/2$ y $1−\alpha/2$ de los datos simulados en cada horizonte de predicción.
La principal ventaja de esta estrategia es que solo requiere de un único modelo para estimar cualquier intervalo. El inconveniente es la necesidad de ejecutar cientos o miles de iteraciones de bootstrapping lo cual resulta muy costoso desde el punto de vista computacional y no siempre es posible.
Librerías¶
# Preprocesado de datos
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset
# Plots
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
pio.renderers.default = 'notebook'
poff.init_notebook_mode(connected=True)
plt.style.use('seaborn-v0_8-darkgrid')
from skforecast.plot import plot_residuals
from skforecast.plot import plot_prediction_distribution
from pprint import pprint
# Modelado y Forecasting
# ==============================================================================
import skforecast
import sklearn
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.preprocessing import OneHotEncoder
from sklearn.compose import ColumnTransformer
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregDirect import ForecasterAutoregDirect
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from sklearn.metrics import mean_pinball_loss
from scipy.stats import norm
# Configuración
# ==============================================================================
import warnings
warnings.filterwarnings('once')
print(f"Version skforecast: {skforecast.__version__}")
print(f"Version sklearn: {sklearn.__version__}")
print(f"Version lightgbm: {lightgbm.__version__}")
Datos¶
# Descarga de datos
# ==============================================================================
data = fetch_dataset(name='bike_sharing_extended_features')
data.head(2)
# One hot encoding de las variables categóricas
# ==============================================================================
encoder = ColumnTransformer(
[('one_hot_encoder', OneHotEncoder(sparse_output=False), ['weather'])],
remainder='passthrough',
verbose_feature_names_out=False
).set_output(transform="pandas")
data = encoder.fit_transform(data)
# Selección de las variables exógenas
# ==============================================================================
exog_features = [
'weather_clear', 'weather_mist', 'weather_rain', 'month_sin', 'month_cos',
'week_of_year_sin', 'week_of_year_cos', 'week_day_sin', 'week_day_cos',
'hour_day_sin', 'hour_day_cos', 'sunrise_hour_sin', 'sunrise_hour_cos',
'sunset_hour_sin', 'sunset_hour_cos', 'temp', 'holiday'
]
data = data[['users'] + exog_features]
Para facilitar el entrenamiento de los modelos, la búsqueda de hiperparámetros óptimos y la evaluación de su capacidad predictiva, los datos se dividen en tres conjuntos separados: entrenamiento, validación y prueba.
# Partición de datos en entrenamiento-validación-test
# ==============================================================================
data = data.loc['2011-05-30 23:59:00':, :]
end_train = '2012-08-30 23:59:00'
end_validation = '2012-11-15 23:59:00'
data_train = data.loc[: end_train, :]
data_val = data.loc[end_train:end_validation, :]
data_test = data.loc[end_validation:, :]
print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Dates validacion : {data_val.index.min()} --- {data_val.index.max()} (n={len(data_val)})")
print(f"Dates test : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
Exploración gráfica¶
La exploración gráfica de series temporales puede ser una forma eficaz de identificar tendencias, patrones y variaciones estacionales. Esto, a su vez, ayuda a orientar la selección del modelo de forecasting más adecuado.
# Gráfico de la serie temporal
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=data_train.index, y=data_train['users'], mode='lines', name='Train'))
fig.add_trace(go.Scatter(x=data_val.index, y=data_val['users'], mode='lines', name='Validation'))
fig.add_trace(go.Scatter(x=data_test.index, y=data_test['users'], mode='lines', name='Test'))
fig.update_layout(
title = 'Número de usuarios',
xaxis_title="Fecha",
yaxis_title="Usuarios",
legend_title="Partición:",
width=900,
height=400,
margin=dict(l=20, r=20, t=35, b=20),
legend=dict(
orientation="h",
yanchor="top",
y=1,
xanchor="left",
x=0.001
)
)
fig.show()
Los gráficos de autocorrelación muestran la correlación entre una serie temporal y sus valores pasados. Son una herramienta útil para identificar el orden de un modelo autorregresivo, es decir, los valores pasados (lags) que se deben incluir en el modelo.
La función de autocorrelación (ACF) mide la correlación entre una serie temporal y sus valores pasados. La función de autocorrelación parcial (PACF) mide la correlación entre una serie temporal y sus valores pasados, pero solo después de eliminar las variaciones explicadas por los valores pasados intermedios.
# Gráfico de autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2))
plot_acf(data.users, ax=ax, lags=24*3)
plt.show()
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2))
plot_pacf(data.users, ax=ax, lags=24*3)
plt.show()
Los resultados del estudio de autocorrelación indican una correlación significativa entre el número de usuarios en las horas anteriores, así como en los días previos. Esto significa que conocer del número de usuarios durante periodos específicos del pasado proporciona información útil para predecir el número de usuarios en el futuro.
Residuos in-sample¶
Se entrena un recursive-multi-step forecaster y se optimiza sus hiperparámetros. Luego, se estiman los intervalos de predicción utilizando bootstrapping.
# Create forecaster and hyperparameters search
# ==============================================================================
# Forecaster
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 7 # Place holder that will be replaced during the search
)
# Lags used as predictors
lags_grid = [24, 48, (1, 2, 3, 23, 24, 25, 47, 48, 49, 71, 72, 73, 364*24, 365*24)]
# Regressor hyperparameters search space
def search_space(trial):
search_space = {
'lags' : trial.suggest_categorical('lags', lags_grid),
'n_estimators' : trial.suggest_int('n_estimators', 200, 800, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 8, step=1),
'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 25, 500),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.5),
'feature_fraction': trial.suggest_float('feature_fraction', 0.5, 0.8, step=0.1),
'max_bin' : trial.suggest_int('max_bin', 50, 100, step=25),
'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1, step=0.1)
}
return search_space
results_search, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = data.loc[:end_validation, 'users'],
exog = data.loc[:end_validation, exog_features],
steps = 24,
metric = 'mean_absolute_error',
search_space = search_space,
initial_train_size = len(data[:end_train]),
refit = False,
n_trials = 20, # Increase for more exhaustive search
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
best_params = results_search['params'].iloc[0]
best_lags = results_search['lags'].iloc[0]
Una vez encontrados los mejores hiperparámetros, se utiliza la función backtesting_forecaster() para generar los intervalos de predicción de todo el conjunto de test.
El argumento
intervalse emplea para especificar la probabilidad de cobertura deseada de los intervalos de predicción. En este caso,intervalse establece en[10, 90], lo que significa que los intervalos de predicción se calculan para los percentiles 10 y 90, lo que da como resultado una cobertura teórica del 80%.El argumento
n_bootse utiliza para especificar el número de muestras bootstrap que se utilizarán para estimar los intervalos de predicción. Cuanto mayor sea el número de muestras, más precisos serán los intervalos de predicción, pero mayor el tiempo necesario.
Por defecto, los intervalos se calculan utilizando los residuos in-sample (residuos del conjunto de entrenamiento). Esto puede dar lugar a intervalos demasiado estrechos (demasiado optimistas).
# Backtesting con intervalos en los datos de test usando residuos in-sample
# ==============================================================================
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data['users'],
exog = data[exog_features],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(data.loc[:end_validation]),
refit = False,
interval = [10, 90], # 80% intervalos
n_boot = 250,
in_sample_residuals = True, # Usar residuos in-sample
binned_residuals = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
predictions.head(5)
# Function to plot predicted intervals
# ======================================================================================
def plot_predicted_intervals(
predictions: pd.DataFrame,
y_true: pd.DataFrame,
target_variable: str,
initial_x_zoom: list=None,
title: str=None,
xaxis_title: str=None,
yaxis_title: str=None,
):
"""
Plot predicted intervals vs real values
Parameters
----------
predictions : pandas DataFrame
Predicted values and intervals.
y_true : pandas DataFrame
Real values of target variable.
target_variable : str
Name of target variable.
initial_x_zoom : list, default `None`
Initial zoom of x-axis, by default None.
title : str, default `None`
Title of the plot, by default None.
xaxis_title : str, default `None`
Title of x-axis, by default None.
yaxis_title : str, default `None`
Title of y-axis, by default None.
"""
fig = go.Figure([
go.Scatter(name='Prediction', x=predictions.index, y=predictions['pred'], mode='lines'),
go.Scatter(name='Real value', x=y_true.index, y=y_true[target_variable], mode='lines'),
go.Scatter(
name='Upper Bound', x=predictions.index, y=predictions['upper_bound'],
mode='lines', marker=dict(color="#444"), line=dict(width=0), showlegend=False
),
go.Scatter(
name='Lower Bound', x=predictions.index, y=predictions['lower_bound'],
marker=dict(color="#444"), line=dict(width=0), mode='lines',
fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
)
])
fig.update_layout(
title=title,
xaxis_title=xaxis_title,
yaxis_title=yaxis_title,
width=900,
height=400,
margin=dict(l=20, r=20, t=35, b=20),
hovermode="x",
xaxis=dict(range=initial_x_zoom),
legend=dict(
orientation="h",
yanchor="top",
y=1.1,
xanchor="left",
x=0.001
)
)
fig.show()
def empirical_coverage(y, lower_bound, upper_bound):
"""
Calculate coverage of a given interval
"""
return np.mean(np.logical_and(y >= lower_bound, y <= upper_bound))
# Gráfico intervalos
# ==============================================================================
plot_predicted_intervals(
predictions=predictions,
y_true=data_test,
target_variable="users",
initial_x_zoom = ['2012-12-01', '2012-12-25'],
title="Valor real vs predicciones en el conjunto de test",
xaxis_title="Date time",
yaxis_title="users",
)
# Cobertura del intervalo en los datos de test
# ==============================================================================
coverage = empirical_coverage(
y = data.loc[end_validation:, 'users'],
lower_bound = predictions["lower_bound"],
upper_bound = predictions["upper_bound"]
)
print(f"Cobertura del intervalo: {round(100*coverage, 2)} %")
# Area del intervalo
# ==============================================================================
area = (predictions["upper_bound"] - predictions["lower_bound"]).sum()
print(f"Area del intervalo: {round(area, 2)}")
Los intervalos de predicción presentan un exceso de confianza, ya que tienden a ser demasiado estrechos, lo que da lugar a una cobertura real inferior a la cobertura teóritca (80%). Esto se debe a la tendencia de los residuos de entrenamiento (in-sample) a sobrestimar la capacidad de predicción del modelo.
Residuos Out-sample (no condicionados a los valores predichos)¶
El método set_out_sample_residuals() se utiliza para almacenar residuos fuera de muestra (out-sample) calculados con un conjunto de validación a través de backtesting. Una vez añadidos los nuevos residuos al foreacster, se tiene que indicar in_sample_residuals en False para que el forecaster los utilice.
# Backtesting con datos de validación para obtener residuos out-sample
# ==============================================================================
_, predictions_val = backtesting_forecaster(
forecaster = forecaster,
y = data.loc[:end_validation, 'users'],
exog = data.loc[:end_validation, exog_features],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(data.loc[:end_train]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
residuals = data.loc[predictions_val.index, 'users'] - predictions_val['pred']
residuals = residuals.dropna()
# Distribución residuos out-sample
# ==============================================================================
print(pd.Series(np.where(residuals < 0, 'negarivo', 'positivo')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(residuals=residuals, figsize=(7, 4))
# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(residuals=residuals)
# Backtesting con intervalos en los datos de test usando residuos out-sample
# ==============================================================================
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data['users'],
exog = data[exog_features],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(data.loc[:end_validation]),
refit = False,
interval = [10, 90], # 80% intervalos
n_boot = 250,
in_sample_residuals = False, # Usar residuos out-sample
binned_residuals = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
predictions.head(5)
# Gráfico intervalos
# ==============================================================================
plot_predicted_intervals(
predictions=predictions,
y_true=data_test,
target_variable="users",
initial_x_zoom = ['2012-12-01', '2012-12-25'],
title="Valor real vs predicciones en el conjunto de test",
xaxis_title="Date time",
yaxis_title="users",
)
# Cobertura del intervalo en los datos de test
# ==============================================================================
coverage = empirical_coverage(
y = data.loc[end_validation:, 'users'],
lower_bound = predictions["lower_bound"],
upper_bound = predictions["upper_bound"]
)
print(f"Cobertura del intervalo: {round(100*coverage, 2)} %")
# Area del intervalo
# ==============================================================================
area = (predictions["upper_bound"] - predictions["lower_bound"]).sum()
print(f"Area del intervalo: {round(area, 2)}")
Los intervalos de predicción obtenidos utilizando los residuos out-sample muestran una mayor amplitud en comparación con los generados utilizando los residuos in-sample, alcanzando una cobertura empírica más próxima a la cobertura nominal. Sin embargo, al examinar el gráfico, es fácil ver que los intervalos son especialmente amplios cuando los valores predichos son bajos, lo que indica que el modelo no es capaz de localizar correctamente la incertidumbre de sus predicciones.
Residuos Out-sample (condicionados a los valores predichos)¶
El proceso de bootstrapping asume que los residuos se distribuyen de forma independiente por lo que pueden utilizarse sin tener en cuenta del valor predicho. En realidad, esto rara vez es cierto; en la mayoría de los casos, la magnitud de los residuos está correlacionado con la magnitud del valor predicho. En este caso, por ejemplo, difícilmente cabría esperar que el error fuera el mismo cuando el número previsto de usuarios es cercano a cero que cuando es de varios cientos.
Para tener en cuenta esta dependencia, skforecast permite distribuir los residuos en K intervalos, en los que cada intervalo está asociado a un rango de valores predichos. Con esta estrategia, el proceso de bootstrapping muestrea los residuos de diferentes intervalos en función del valor predicho, lo que puede mejorar la cobertura del intervalo y ajustar su anchura cuando sea necesario, permitiendo que el modelo distribuya mejor la incertidumbre de sus predicciones.
Para que el forecaster pueda agrupar los residuos out-sample, los valores predichos se pasan al método set_out_sample_residuals() junto con los residuos. Internamente, skforecast utiliza un sklearn.preprocessing.KBinsDiscretizer con los parámetros: n_bins=15, encode=ordinal, strategy=quantile,
subsample=10000, random_state=789654, dtype=np.float64. El proceso de binning puede ajustarse utilizando el argumento binner_kwargs del objeto forecaster.
# Crear forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1, **best_params),
lags = best_lags,
binner_kwargs = {'n_bins': 10}
)
forecaster.fit(
y = data.loc[:end_validation, 'users'],
exog = data.loc[:end_validation, exog_features]
)
Durante el proceso de entrenamiento, el forecaster utiliza las predicciones in sample para definir los intervalos en los que se almacenan los residuos, en función del valor prediccho con el que están relacionados. Aunque no se utilizan en este ejemplo, los residuos in-sample también se dividen en intervalos y se almacenan en el atributo in_sample_residuals_by_bin.
# Intervalos de los bins
# ==============================================================================
pprint(forecaster.binner_intervals)
A continuación, los residuos out-sample se almacenan en el forecaster. Dado que también se proporcionan los valores predichos, se clasifican según los intervalos aprendidos durante el entrenamiento. Para evitar un uso excesivo de memoria, se almacena un máximo de 200 residuos por partición (bin).
# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(residuals=residuals, y_pred=predictions_val['pred'])
# Número de residuos por bin
# ==============================================================================
for k, v in forecaster.out_sample_residuals_by_bin.items():
print(f" Bin {k}: n={len(v)}")
# Distribución de los residuos por bin
# ==============================================================================
out_sample_residuals_by_bin_df = pd.DataFrame(
dict([(k, pd.Series(v)) for k,v in forecaster.out_sample_residuals_by_bin.items()])
)
fig, ax = plt.subplots(figsize=(6, 3))
out_sample_residuals_by_bin_df.boxplot(ax=ax)
ax.set_title("Distribución de los residuos por bin")
ax.set_xlabel("Bin")
ax.set_ylabel("Resoduals");
Los gráficos de caja muestran cómo difieren la dispersión y la magnitud de los residuos en función del valor predicho. Por ejemplo, para el bin 0, los residuos nunca superan un valor absoluto de 100, mientras que para el bin 9, para valores predichos en el intervalo a menudo sí lo hacen.
Por último, los intervalos de predicción para los datos de test se estiman mediante el proceso de backtesting, con los residuos out-sample condicionados a los valores predichos.
# Backtesting con intervalos en los datos de test usando residuos out-sample
# ==============================================================================
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data['users'],
exog = data[exog_features],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(data.loc[:end_validation]),
refit = False,
interval = [10, 90], # 80% intervalos
n_boot = 250,
in_sample_residuals = False, # Usar residuos out-sample
binned_residuals = True, # Usar residuos por bin
n_jobs = 'auto',
verbose = False,
show_progress = True
)
predictions.head(5)
# Gráfico intervalos
# ==============================================================================
plot_predicted_intervals(
predictions=predictions,
y_true=data_test,
target_variable="users",
initial_x_zoom = ['2012-12-01', '2012-12-25'],
title="Valor real vs predicciones en el conjunto de test",
xaxis_title="Date time",
yaxis_title="users",
)
# Cobertura del intervalo en los datos de test
# ==============================================================================
coverage = empirical_coverage(
y = data.loc[end_validation:, 'users'],
lower_bound = predictions["lower_bound"],
upper_bound = predictions["upper_bound"]
)
print(f"Cobertura del intervalo: {round(100*coverage, 2)} %")
# Area del intervalo
# ==============================================================================
area = (predictions["upper_bound"] - predictions["lower_bound"]).sum()
print(f"Area del intervalo: {round(area, 2)}")
Cuando se utilizan residuos out-sample condicionados al valor predicho, el intervalo tiene una cobertura próxima al valor esperado (80%) a la vez que se reduce su amplitud. El modelo es capaz de distribuir mejor la incertidumbre en sus predicciones.
Predicción: bootstraping, cuantiles y distribución¶
En las secciones anteriores se ha mostrado el uso del proceso de backtesting para estimar el intervalo de predicción a lo largo de un periodo de tiempo determinado. El objetivo es imitar el comportamiento del modelo en producción ejecutando predicciones a intervalos regulares, actualizando incrementalmente los datos de entrada.
También, es posible ejecutar una única predicción que estime N strps por delante sin pasar por todo el proceso de backtesting. En estos casos, skforecast ofrece cuatro métodos diferentes: predict_bootstrapping, predict_interval, predict_quantile y predict_distribution.
Predict Bootstraping
El método predict_bootstrapping realiza las iteraciones de bootstrapping para generar la muestra de posibles predicciones alternativas. Estos son los valores subyacentes utilizados para calcular los intervalos, cuantiles y distribuciones.
# Entrenar forecaster
# ==============================================================================
forecaster.fit(
y = data.loc[:end_validation, 'users'],
exog = data.loc[:end_validation, exog_features]
)
# Predecir 10 secuencias diferentes de 7 steps cada una usando bootstrapping
# ==============================================================================
boot_predictions = forecaster.predict_bootstrapping(
exog = data_test[exog_features],
steps = 7,
n_boot = 25
)
boot_predictions
Un gráfico de ridge es una forma útil de visualizar la incertidumbre asociada a las predicciones de un modelo de forecasting. Este gráfico muestra la distribución de desnidad para cada step utilizando las predicciones obtenidas por bootstrapping.
# Ridge plot de las distribuciones de bootstrapping
# ==============================================================================
_ = plot_prediction_distribution(boot_predictions, figsize=(7, 4))
Predict Interval
En la mayoría de los casos, el usuario está interesado en un intervalo específico en lugar de toda la matriz de predicciones obtenida por bootstrapping. Para cubrir esta necesidad, skforecast proporciona el método predict_interval. Este método utiliza internamente predict_bootstrapping y estima los cuantiles superior e inferior para cada step, proporcionando así los intervalos de predicción deseados.
# Predicción de intervalos para los siguiente 7 steps
# ==============================================================================
predictions = forecaster.predict_interval(
exog = data_test[exog_features],
steps = 7,
interval = [10, 90],
n_boot = 150
)
predictions
Predict Quantile
Este método funciona de forma similar a predict_interval, con la característica añadida de que permite definir una lista específica de cuantiles. Es importante recordar que estos cuantiles deben estar dentro del intervalo [0, 1].
# Predicción de cuantiles 5, 25, 75 and 95 para los siguientes 7 steps
# ==============================================================================
predictions = forecaster.predict_quantiles(
exog = data_test[exog_features],
steps = 7,
n_boot = 150,
quantiles = [0.05, 0.25, 0.75, 0.95],
)
predictions
Predict Distribution
Las estimaciones de intervalos mostradas en los apartados anteriores no hacen ninguna asunción sobre una distribución concreta. El método predict_dist de skforecast permite ajustar una distribución paramétrica a las simulaciones obtenidas con predict_bootstrapping. Esto es útil cuando hay razones para creer que los errores de predicción siguen una distribución particular, como la distribución normal o la distribución t-student. El método predict_dist permite al usuario especificar cualquier distribución continua del módulo scipy.stats.
# Predicción de los parámetros de una distribución normal en cada uno de los siguentes 7 steps
# ==============================================================================
predictions = forecaster.predict_dist(
exog = data_test[exog_features],
steps = 7,
n_boot = 150,
distribution = norm
)
predictions
Intervalos de predicción utilizando modelos de regresión cuantílica¶
A diferencia de los modelos de regresión más comunes, que pretende estimar la media de la variable respuesta dados ciertos valores de las variables predictoras, la regresión cuantílica tiene como objetivo estimar los cuantiles condicionales de la variable respuesta. Para una función de distribución continua, el cuantil $\alpha$ $Q_{\alpha}(x)$ se define como el valor tal que la probabilidad de que $Y$ sea menor que $Q_{\alpha}(x)$ es, para un determinado $X=x$, igual a $\alpha$. Por ejemplo, el 36% de los valores de la población son inferiores al cuantil $Q=0,36$. El cuantil más conocido es el cuantil 50%, más comúnmente conocido como mediana.
Al combinar las predicciones de dos modelos de regresión cuantílica, es posible construir un intervalo donde, cada modelo, estima uno de los límites del intervalo. Por ejemplo, los modelos obtenidos para $Q = 0.1$ y $Q = 0.9$ generan un intervalo de predicción del 80% (90% - 10% = 80%).
Son varios los algoritmos de machine learning capaces de modelar cuantiles. Algunos de ellos son:
Así como el error cuadrático se utiliza como función de coste para entrenar modelos que predicen el valor medio, se necesita una función de coste específica para entrenar modelos que predicen cuantiles. La función utilizada con más frecuencia para la regresión de cuantiles se conoce como pinball:
$$\text{pinball}(y, \hat{y}) = \frac{1}{n_{\text{muestras}}} \sum_{i=0}^{n_{\text{muestras}}- 1} \alpha \max(y_i - \hat{y}_i, 0) + (1 - \alpha) \max(\hat{y}_i - y_i, 0)$$donde $\alpha$ es el cuantil objetivo, $y$ el valor real e $\hat{y}$ la predicción del cuantil.
Se puede observar que el coste difiere según el cuantil evaluado. Cuanto mayor sea el cuantil, más se penalizan las subestimaciones y menos las sobreestimaciones. Al igual que con MSE y MAE, el objetivo es minimizar sus valores (a menor coste, mejor).
Dos desventajas de la regresión por cuantiles en comparación con el método de bootstrapping son: que cada cuantil necesita su regresor, y que la regresión por cuantiles no está disponible para todos los tipos de modelos de regresión. Sin embargo, una vez entrenados los modelos, la inferencia es mucho más rápida ya que no se necesita un proceso iterativo.
Este tipo de intervalos de predicción se pueden estimar fácilmente utilizando ForecasterAutoregDirect.
⚠ Warning
Los forecaters de tipoForecasterAutoregDirect son más lentos que los ForecasterAutoregRecursive porque requieren el entrenamiento de un modelo por *step*. Para limitar el tiempo necesario para ejecutar los siguientes ejemplos, los datos se agregan a frecuencia diaria y sólo se predicen 7 los siguientes 7 dias (una semana).
Datos¶
# Descarga de datos
# ==============================================================================
data = fetch_dataset(name='bike_sharing_extended_features', verbose=False)
# Agragar datos a frecuencia horaria
# ==============================================================================
data = (
data
.resample(rule="D", closed="left", label="right")
.agg({"users": "sum"})
)
# Particion train-validation-test
# ==============================================================================
end_train = '2012-05-31 23:59:00'
end_validation = '2012-09-15 23:59:00'
data_train = data.loc[: end_train, :]
data_val = data.loc[end_train:end_validation, :]
data_test = data.loc[end_validation:, :]
print(f"Fechas train : {data_train.index.min()} --- {data_train.index.max()} (n={len(data_train)})")
print(f"Fechas validacion : {data_val.index.min()} --- {data_val.index.max()} (n={len(data_val)})")
print(f"Fechas test : {data_test.index.min()} --- {data_test.index.max()} (n={len(data_test)})")
Modelos de regresión por cuantiles¶
Se estima un intervalo de predicción del 80% para los siguientes 7 días. En este ejemplo, se entrena un modelo de gradient boosting LightGBM, sin embargo, el lector puede utilizar cualquier otro modelo simplemente reemplazando la definición del regresor.
# Crear forecasters: uno para cada límite del intervalo
# ==============================================================================
# Los forecasters obtenidos para alpha=0.1 y alpha=0.9 producen un intervalo de
# confianza del 80% (90% - 10% = 80%).
# Forecaster para cuantil 10%
forecaster_q10 = ForecasterAutoregDirect(
regressor = LGBMRegressor(
objective = 'quantile',
metric = 'quantile',
alpha = 0.1,
random_state = 15926,
verbose = -1
),
lags = 7,
steps = 7
)
# Forecaster para cuantil 90%
forecaster_q90 = ForecasterAutoregDirect(
regressor = LGBMRegressor(
objective = 'quantile',
metric = 'quantile',
alpha = 0.9,
random_state = 15926,
verbose = -1
),
lags = 7,
steps = 7
)
Al validar un modelo de regresión por cuantiles, se debe proporcionar una métrica personalizada según el cuantil que se esté estimando.
# Función de coste para cada cuantil (pinball_loss)
# ==============================================================================
def mean_pinball_loss_q10(y_true, y_pred):
"""
Pinball loss for quantile 10.
"""
return mean_pinball_loss(y_true, y_pred, alpha=0.1)
def mean_pinball_loss_q90(y_true, y_pred):
"""
Pinball loss for quantile 90.
"""
return mean_pinball_loss(y_true, y_pred, alpha=0.9)
# Grid search para los hyper-parámetros y lags de cada quantile-forecaster
# ==============================================================================
def search_space(trial):
search_space = {
'n_estimators' : trial.suggest_int('n_estimators', 100, 500, step=50),
'max_depth' : trial.suggest_int('max_depth', 3, 10, step=1),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.1)
}
return search_space
results_grid_q10 = bayesian_search_forecaster(
forecaster = forecaster_q10,
y = data.loc[:end_validation, 'users'],
steps = 7,
metric = mean_pinball_loss_q10,
search_space = search_space,
initial_train_size = len(data[:end_train]),
refit = False,
n_trials = 10, # Increase for more exhaustive search
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
results_grid_q90 = bayesian_search_forecaster(
forecaster = forecaster_q90,
y = data.loc[:end_validation, 'users'],
steps = 7,
metric = mean_pinball_loss_q90,
search_space = search_space,
initial_train_size = len(data[:end_train]),
refit = False,
n_trials = 10, # Increase for more exhaustive search
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
Una vez que se han encontrado los mejores hiperparámetros para cada forecaster, se aplica un proceso de backtesting utilizando los datos de test.
# Backtesting con datos de test
# ==============================================================================
metric_q10, predictions_q10 = backtesting_forecaster(
forecaster = forecaster_q10,
y = data['users'],
steps = 7,
metric = mean_pinball_loss_q10,
initial_train_size = len(data.loc[:end_validation]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
metric_q90, predictions_q90 = backtesting_forecaster(
forecaster = forecaster_q90,
y = data['users'],
steps = 7,
metric = mean_pinball_loss_q90,
initial_train_size = len(data.loc[:end_validation]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Gráfico
# ==============================================================================
fig = go.Figure([
go.Scatter(name='Real value', x=data_test.index, y=data_test['users'], mode='lines'),
go.Scatter(
name='Upper Bound', x=predictions_q90.index, y=predictions_q90['pred'],
mode='lines', marker=dict(color="#444"), line=dict(width=0), showlegend=False
),
go.Scatter(
name='Lower Bound', x=predictions_q10.index, y=predictions_q10['pred'],
marker=dict(color="#444"), line=dict(width=0), mode='lines',
fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
)
])
fig.update_layout(
title="Valores reales vs predicciones en los datos de test",
xaxis_title="Date time",
yaxis_title="users",
width=900,
height=400,
margin=dict(l=20, r=20, t=35, b=20),
hovermode="x",
legend=dict(
orientation="h",
yanchor="top",
y=1.1,
xanchor="left",
x=0.001
)
)
fig.show()
# Cobertura del intervalo en los datos de test
# ==============================================================================
coverage = empirical_coverage(
y = data.loc[end_validation:, 'users'],
lower_bound = predictions_q10["pred"],
upper_bound = predictions_q90["pred"]
)
print(f"Cobertura del intervalo: {round(100*coverage, 2)} %")
# Area del intervalo
# ==============================================================================
area = (predictions_q90["pred"] - predictions_q10["pred"]).sum()
print(f"Area del intervalo: {round(area, 2)}")
En este caso de uso, la estrategia de previsión cuantílica no logra una cobertura empírica cercana a la esperada (80%).
Información de sesión¶
import session_info
session_info.show(html=False)
Bibliografía¶
Instrucciones para citar¶
¿Cómo citar este documento?
Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!
Forecasting probabilistico con machine learning 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://www.cienciadedatos.net/documentos/py42-foreasting-probabilistico.html
¿Cómo citar skforecast?
Si utilizas skforecast en tu investigación o publicación, te lo agradeceríamos mucho que lo cites. ¡Muchas gracias!
Zenodo:
Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2024). skforecast (v0.12.0). Zenodo. https://doi.org/10.5281/zenodo.8382788
APA:
Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.12.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.12.0}, month = {4}, 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
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este documento creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.