Gradient boosting probabilístico (NGBoost) con python
Introducción¶
La mayoría de métodos de machine learning aplicados a problemas de regresión modelan el comportamiento de una variable respuesta $y$ en función de uno o varios predictores $X$, y generan predicciones de tipo "estimación puntual". Estás predicciones se corresponden el valor esperado de $y$ dado un determinado valor de los predictores ($\mathbb{E}[y|X]$).
La principal limitación de este tipo de modelos es que no ofrecen ninguna información sobre la incertidumbre asociada al valor predicho. Tampoco permiten responder a preguntas como:
¿Cuál es la probabilidad de que, dado un valor de $X$, la variable respuesta $y$ sea menor o mayor que un determinado valor?
¿Cuál es la probabilidad de que, dado un valor de $X$, la variable respuesta $y$ esté comprendida dentro de un intervalo?
Para responder a estas preguntas, se necesitan métodos de regresión probabilística capaces de modelar la distribución de probabilidad condicional de la variable respuesta en función de los predictores ($P(y|X)$).
Véase el siguiente ejemplo simulado (y muy simplificado) sobre la evolución del consumo eléctrico de todas las casas de una ciudad en función de la hora del día.
# Simulación distribución no uniforme en el rango X
# ==============================================================================
import numpy as np
from scipy.stats import norm
from matplotlib import style
import matplotlib.pyplot as plt
np.random.seed(seed=1234)
n = 3000
x = np.linspace(start=0, stop= 24, num=n)
y = np.random.normal(
loc = 15,
scale = 1 + 1.5*((4.8 < x) & (x < 7.2)) + 4*((7.2 < x) & (x < 12)) \
+ 1.5*((12 < x) & (x < 14.4)) + 2*(x > 16.8)
)
# Gráfico
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))
ax.hlines(y=15, xmin=0, xmax=24, color='#FC4E07', label='media')
ax.scatter(x, y, alpha = 0.1, c = 'k')
ax.set_xticks(range(0,25))
ax.set_title('Evolución del consumo eléctrico a lo largo del día', fontdict={'fontsize':13})
ax.set_xlabel('Hora del día')
ax.set_ylabel('Consumo eléctrico (MWh)')
plt.legend();
La media del consumo eléctrico es la misma durante todo el día, $\overline{consumo}$ = 15 Mwh, sin embargo, su dispersión no es constante (heterocedasticidad). Véase el resultado de predecir el consumo medio en función de la hora del día con un modelo Gradient Boosting.
El valor predicho es muy próximo a la media real, es decir, el modelo es bueno prediciendo el consumo medio esperado. Ahora, imagínese que la compañía encargada de suministrar la electricidad debe de ser capaz de provisionar, en un momento dado, con hasta un 50% de electricidad extra respecto al promedio. Esto significa un máximo de 22.5 Mwh. Estar preparado para suministrar este extra de energía implica gastos de personal y maquinaría, por lo que la compañía se pregunta si es necesario estar preparado para producir tal cantidad durante todo el día, o si, por lo contrario, podría evitarse durante algunas horas, ahorrando así gastos.
Un modelo que predice únicamente el promedio no permite responder a esta pregunta ya que, tanto para las 2h de la mañana como para las 8h, el consumo promedio predicho es en torno a 15 Mwh. Sin embargo, la probabilidad de que se alcancen consumos de 22.5 Mwh a las 2h es prácticamente nula mientras que esto ocurra a las 10h sí es razonable.
Una forma de poder predecir la probabilidad de que el consumo eléctrico exceda un determinado valor, a una hora determinada, es modelar su distribución de probabilidad en función de la hora del día.
Natural Gradient Boosting (NGBoost)¶
Natural Gradient Boosting (NGBoost) es una generalización de gradient boosting que capaz de estimar los parámetros $\theta$ de una distribución de probabilidad condicional, permitiendo así generar predicciones probabilísticas. Sin entrar en detalles matemáticos, las dos principales diferencias respecto a gradient boosting son:
Utiliza el descenso de gradiente natural en lugar del descenso de gradiente para el ajuste del modelo.
Estima los parámetros de una distribución de probabilidad paramétrica $P_{\theta}(y|X)$ en lugar de una estimación puntual $\mathbb{E}[y|X]$.
Por ejemplo, si la distribución elegida es una normal, NGBoost modela los valores de la media $\mu$ y desviación típica $\sigma$ en función de los predictores $X$. De esta forma, no es necesario asumir la condición de homocedasticidad (misma varianza para todo el rango de $X$).
Natural Gradient Boosting con Python
La librería ngboost desarrollada por el equipo Stanford Machine Learning Group implementa este tipo de modelos en python, con una API similar y compatible con la de scikit learn.
En términos generales, para entrenar un modelo ngboost es necesario definir 3 elementos:
Un tipo de base learner, por defecto se emplea árboles de decisión de scikit-learn.
Una distribución paramétrica de probabilidad.
Una regla de scoring que permita dirigir el ajuste.
El resto de pasos son iguales a los seguidos en un modelo gradient boosting de scikitlearn, incluida la posibilidad de optimizar los hiperparámetros mediante GridSearchCV y early_stopping_rounds.
La librería ngboost es reciente y seguro que irá evolucionando e incluyendo más funcionalidades. Puede encontrarse información detallada en NGBoost user guide. Por el momento, las distribuciones disponibles y reglas de scoring son:
Ejemplo 1¶
Una empresa suministradora de energía eléctrica dispone del consumo horario de todas las casas de una ciudad. Para distribuir sus recursos de la forma más óptima, quiere ser capaz generar un modelo probabilístico que le permita predecir, con un intervalo de confianza del 90%, el consumo promedio para cada hora del día. También quieren poder calcular la probabilidad a lo largo del día de que el consumo exceda un determinado valor.
# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style
import seaborn as sns
# Modelado
# ==============================================================================
from ngboost import NGBRegressor
from ngboost import scores
from ngboost import distns
from sklearn.tree import DecisionTreeRegressor
from scipy import stats
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
Datos¶
Se simulan datos a partir de una distribución normal con una varianza distinta dependiendo de la hora del día.
# Datos simulados
# ==============================================================================
# Simulación ligeramente modificada del ejemplo publicado en
# XGBoostLSS – An extension of XGBoost to probabilistic forecasting Alexander März
np.random.seed(seed=1234)
n = 3000
# Distribución normal con varianza cambiante
x = np.linspace(start=0, stop= 24, num=n)
y = np.random.normal(
loc = 15,
scale = 1 + 1.5*((4.8 < x) & (x < 7.2)) + 4*((7.2 < x) & (x < 12)) \
+ 1.5*((12 < x) & (x < 14.4)) + 2*(x > 16.8)
)
Dado que los datos se han simulado empleando distribuciones normales, se conoce el valor de los cuantiles teóricos para cada hora.
# Cálculo del cuantil 0.1 y 0.9 para cada posición de x simulada
# ==============================================================================
cuantil_10 = stats.norm.ppf(
q = 0.1,
loc = 15,
scale = 1 + 1.5*((4.8 < x) & (x < 7.2)) + 4*((7.2 < x) & (x < 12)) \
+ 1.5*((12 < x) & (x < 14.4)) + 2*(x > 16.8)
)
cuantil_90 = stats.norm.ppf(
q = 0.9,
loc = 15,
scale = 1 + 1.5*((4.8 < x) & (x < 7.2)) + 4*((7.2 < x) & (x < 12)) \
+ 1.5*((12 < x) & (x < 14.4)) + 2*(x > 16.8)
)
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))
ax.scatter(x, y, alpha = 0.2, c = "#333")
ax.plot(x, cuantil_10, c = "black")
ax.plot(x, cuantil_90, c = "black")
ax.fill_between(x, cuantil_10, cuantil_90, alpha=0.2, color='red',
label='Intevalo cuantílico teórico 0.1-0.9')
ax.set_xticks(range(0,25))
ax.set_title('Evolución del consumo eléctrico a lo largo del día',
fontdict={'fontsize':13})
ax.set_xlabel('Hora del día')
ax.set_ylabel('Consumo eléctrico (MWh)')
plt.legend();
Modelo¶
Para calcular los cuantiles, es necesario obtener primero la distribución probabilística de los datos. Con el objetivo de facilitar la demostración de este ejemplo, se asumen dos simplificaciones:
Se asume que la distribución es de tipo normal (se han simulado con esta distribución).
Se utilizan los hiperparámetros por defecto al ajustar el modelo.
En la práctica, es necesario un análisis preliminar que permita determinar estos dos puntos (ver ejemplo 2).
# Ajuste del modelo
# ==============================================================================
modelo = NGBRegressor(
Dist = distns.Normal,
Base = DecisionTreeRegressor(criterion='friedman_mse', max_depth=3),
Score = scores.LogScore,
natural_gradient = True,
n_estimators = 500,
learning_rate = 0.01,
minibatch_frac = 1.0,
col_sample = 1.0,
random_state = 123,
)
modelo.fit(X=x.reshape(-1, 1), Y=y)
Predicciones¶
Los modelos NGBRegressor tienen dos métodos de predicción:
predict(): devuelve una predicción de tipo puntual (point estimated), tal como hacen la mayoría de modelos. Este es el valor esperado de la variable respuesta acorde al valor de los predictores en la observación de test $\mathbb{E}[y|x=x_i]$.pred_dist(): devuelve un objetongboost.distnsque contiene la distribución condicional de la variable respuesta para el valor de los predictores en la observación de test $P(y|x=x_i)$. Dentro de cada uno de estos objetos se encuentran almacenados los parámetros estimados.
# Predicción point estimated
# ==============================================================================
predicciones_point = modelo.predict(X=x.reshape(-1, 1))
predicciones_point
# Predicción distribución condicional
# ==============================================================================
predicciones_dist = modelo.pred_dist(X=x.reshape(-1, 1))
predicciones_dist
# Dentro de cada objeto distns.distn se encuentran los parámetros estimados
predicciones_dist[0].params
# Parámetros estimados para cada hora
# ==============================================================================
df_parametros = pd.DataFrame([dist.params for dist in predicciones_dist])
df_parametros['hora'] = x
df_parametros
Predicción intervalos¶
Una vez conocidos los parámetros de la distribución condicional, se pueden calcular los cuantiles del consumo energético esperado para cada hora.
Un cuantil de orden $\tau$ $(0 < \tau < 1)$ de una distribución es el valor que marca un corte tal que, una proporción $\tau$ de valores de la población, es menor o igual que dicho valor. Por ejemplo, el cuantil de orden 0.36 deja un 36% de valores por debajo y el cuantil de orden 0.50 el 50% (se corresponde con la mediana de la distribución).
Todas las distribuciones de generadas por NGBRegressor son clases de scipy.stats, por lo que disponen de los métodos pdf(), logpdf(), cdf(), ppf() y rvs() con los que calcular la densidad, logaritmo de densidad, probabilidad acumulada, cuantiles, y muestreo de nuevos valores.
Los pasos a seguir para calcular los intervalos cuantílicos son:
Generar un grid de valores que cubra todo el intervalo de valores observado.
Para cada valor del grid, obtener los parámetros predichos por el modelo.
Calcular los cuantiles de la distribución seleccionada con los parámetros predichos.
# Cálculo de los cuantiles teóricos para establecer el intervalo central que
# acumula un 90% de probabilidad empleando los parámetros predichos.
# ==============================================================================
grid_X = np.linspace(start=min(x), stop=max(x), num=1000)
predicciones_dist = modelo.pred_dist(X=grid_X.reshape(-1, 1))
pred_cuantil_10 = [stats.norm.ppf(q=0.1, **dist.params) for dist in predicciones_dist]
pred_cuantil_90 = [stats.norm.ppf(q=0.9, **dist.params) for dist in predicciones_dist]
# Gráfico intervalo cuantílico 10%-90%
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))
ax.scatter(x, y, alpha = 0.2, c = "#333")
ax.plot(x, cuantil_10, c = "black")
ax.plot(x, cuantil_90, c = "black")
ax.fill_between(x, cuantil_10, cuantil_90, alpha=0.2, color='red',
label='Intevalo cuantílico teórico 0.1-0.9')
ax.plot(grid_X, pred_cuantil_10, c = "blue", label='intevalo cuantilico predicho')
ax.plot(grid_X, pred_cuantil_90, c = "blue")
ax.set_xticks(range(0,25))
ax.set_title('Evolución del consumo eléctrico a lo largo del día',
fontdict={'fontsize':13})
ax.set_xlabel('Hora del día')
ax.set_ylabel('Consumo eléctrico (MWh)')
plt.legend();
Cobertura¶
Una de las métricas empleadas para evaluar intervalos es la cobertura (coverage). Su valor se corresponde con el porcentaje de observaciones que caen dentro del intervalo estimado. Idealmente, la cobertura debe de ser igual a la confianza del intervalo, por lo que, en la práctica, cuanto más próximo sea su valor, mejor.
En este ejemplo, se han calculado los cuantiles 0.1 y 0.9, así que, el intervalo, tiene una confianza del 80%. Si el intervalo predicho es correcto, aproximadamente el 80% de las observaciones estarán dentro.
# Predicción intervalos para las observaciones disponibles
# ==============================================================================
predicciones_dist = modelo.pred_dist(X=x.reshape(-1, 1))
pred_cuantil_10 = [stats.norm.ppf(q=0.1, **dist.params) for dist in predicciones_dist]
pred_cuantil_90 = [stats.norm.ppf(q=0.9, **dist.params) for dist in predicciones_dist]
# Cobertura del intervalo predicho
# ==============================================================================
dentro_intervalo = np.where((y >= pred_cuantil_10) & (y <= pred_cuantil_90), True, False)
cobertura = dentro_intervalo.mean()
print(f"Cobertura del intervalo predicho: {100 * cobertura}")
Predicción probabilidad¶
Si por ejemplo, se desea conocer la probabilidad de que a las 8h el consumo supere los 22.5 MWh, primero se predicen los parámetros de la distribución a para $hora=8$ y después, con esos parámetros, se calcula la probabilidad de $consumo≥22.5$ con la función de distribución acumulada.
# Predicción de los parámetros de la distribución cuando x=8
# ==============================================================================
distribucion = modelo.pred_dist(X=np.array([8]).reshape(-1, 1))
distribucion.params
# Probabilidad de que el consumo supere un valor de 22.5
# ==============================================================================
100 * (1 - stats.norm.cdf(x=22.5, **distribucion.params))
Acorde a la distribución predicha por el modelo, existe una probabilidad del 3.24% de que el consumo a las 8h sea igual o superior a 22.5 MWh.
Si este mismo cálculo se realiza para todo el rango horario, se puede conocer la probabilidad de que el consumo exceda un determinado valor a lo largo del día.
# Predicción de probabilidad en todo el rango horario
# ==============================================================================
grid_X = np.linspace(start=0, stop=24, num=500)
predicciones_dist = modelo.pred_dist(X=grid_X.reshape(-1, 1))
predicciones_proba = [(1 - stats.norm.cdf(x=22.5, **dist.params)) for dist in predicciones_dist]
# Gráfico probabilidades
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 4))
ax.plot(grid_X, predicciones_proba)
ax.set_xticks(range(0, 25))
ax.set_title('Probabilidad de que el consumo supere los 22.5 KWh',
fontdict={'fontsize':13})
ax.set_xlabel('Hora del día')
ax.set_ylabel('Probabilidad');