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');
Ejemplo 2¶
El set de datos dbbmi (Fourth Dutch Growth Study, Fredriks et al. (2000a, 2000b), disponible en el paquete de R gamlss.data, contiene información sobre la edad (age) e índice de masa corporal (bmi) de 7294 jóvenes holandeses de entre 0 y 20 años. El objetivo es entrenar un modelo capaz de predecir cuantiles del índice de masa corporal en función de la edad, ya que este es uno de los estándares empleados para detectar casos anómalos que pueden requerir atención médica.
Librerías¶
# 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 sklearn.model_selection import GridSearchCV
from scipy import stats
import inspect
import multiprocessing
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
Datos¶
url = 'https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine-learning-python/master/data/dbbmi.csv'
datos = pd.read_csv(url)
datos.info()
# Gráfico distribución
# ==============================================================================
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(8, 3))
ax.hist(datos.bmi, bins=30, density=True, color='#3182bd', alpha=0.8)
ax.set_title('Distribución valores índice masa corporal')
ax.set_xlabel('bmi')
ax.set_ylabel('densidad');
fig, axs = plt.subplots(ncols=2, nrows=1, figsize=(10, 4.5), sharey=True)
datos[datos.age <= 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[0]
)
axs[0].set_ylim([10, 30])
axs[0].set_title('Edad <= 2.5 años')
datos[datos.age > 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[1]
)
axs[1].set_title('Edad > 2.5 años')
fig.tight_layout()
plt.subplots_adjust(top = 0.85)
fig.suptitle('Distribución índice masa corporal en función de la edad', fontsize = 14);
Esta distribución muestra tres características importantes que hay que tener en cuenta a la hora de modelarla:
La relación entre la edad y el índice de masa corporal no es lineal ni constante. Tiene una relación positiva notable hasta los 0.25 años, después se estabiliza hasta los 10 años y vuelve a ascender notablemente de los 10 a los 20 años.
La varianza es heterogénea (heterocedasticidad), siendo esta menor en edades bajas mayor en edades altas.
La distribución de la variable respuesta no es de tipo normal, muestra asimetría y una cola positiva.
Dadas estas características, se necesita un modelo que:
Sea capaz de aprender relaciones no lineales.
Sea capaz de modelar explícitamente la varianza en función de los predictores, ya que esta no es constante.
Sea capaz de aprender distribuciones asimétricas con una marcada cola positiva.
Selección distribución paramétrica¶
Los modelos NGBoost requieren que se especifique un tipo de distribución paramétrica, por lo que el primer paso es identificar cuál de las tres distribuciones disponibles en NGBoost (normal, lognormal y exponencial) se ajusta mejor a los datos.
Se procede a ajustar cada una de las distribuciones y se evalúan gráficamente superponiéndolas con el histograma. Puede encontrarse información detallada de cómo comparar distribuciones en Ajuste y selección de distribuciones con Python.
def plot_multiple_distribuciones(x, nombre_distribuciones, ax=None):
'''
Esta función ajusta y superpone las curvas de densidad de varias distribuciones
con el histograma de los datos.
Parameters
----------
x : array_like
datos con los que ajustar la distribución.
nombre_distribuciones : list
lista con nombres de distribuciones disponibles en `scipy.stats`.
Returns
-------
resultados: matplotlib.ax
gráfico creado
'''
if ax is None:
fig, ax = plt.subplots(figsize=(7,4))
ax.hist(x=x, density=True, bins=30, color="#3182bd", alpha=0.5)
ax.plot(x, np.full_like(x, -0.01), '|k', markeredgewidth=1)
ax.set_title('Ajuste distribuciones')
ax.set_xlabel('x')
ax.set_ylabel('Densidad de probabilidad')
for nombre in nombre_distribuciones:
distribucion = getattr(stats, nombre)
parametros = distribucion.fit(data=x)
nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
if not p=='x'] + ["loc","scale"]
parametros_dict = dict(zip(nombre_parametros, parametros))
log_likelihood = distribucion.logpdf(x, *parametros).sum()
aic = -2 * log_likelihood + 2 * len(parametros)
bic = -2 * log_likelihood + np.log(x.shape[0]) * len(parametros)
x_hat = np.linspace(min(x), max(x), num=100)
y_hat = distribucion.pdf(x_hat, *parametros)
ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
ax.legend();
return ax
fig, ax = plt.subplots(figsize=(8, 4))
plot_multiple_distribuciones(
x=datos.bmi.to_numpy(),
nombre_distribuciones=['lognorm', 'norm', 'expon'],
ax=ax
);
La representación gráfica muestra claras evidencias de que, entre las distribuciones disponibles, la lognormal es la que mejor se ajusta.
Modelo¶
A diferencia del ejemplo anterior, en este caso, se recurre a validación cruzada para identificar los hiperparámetros óptimos.
X = datos.age.to_numpy()
y = datos.bmi.to_numpy()
# Grid de hiperparámetros
# ==============================================================================
b1 = DecisionTreeRegressor(criterion='friedman_mse', max_depth=2)
b2 = DecisionTreeRegressor(criterion='friedman_mse', max_depth=4)
b3 = DecisionTreeRegressor(criterion='friedman_mse')
param_grid = {
'Base': [b1, b2],
'n_estimators': [100, 500, 1000]
}
# Búsqueda por validación cruzada
# ==============================================================================
grid = GridSearchCV(
estimator = NGBRegressor(
Dist = distns.LogNormal,
Score = scores.LogScore,
verbose = False
),
param_grid = param_grid,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = 5,
verbose = 0
)
# Se asigna el resultado a _ para que no se imprima por pantalla
_ = grid.fit(X = X.reshape(-1, 1), y = y)
# Resultados del grid
# ==============================================================================
resultados = pd.DataFrame(grid.cv_results_)
resultados.filter(regex = '(param.*|mean_t|std_t)')\
.drop(columns = 'params')\
.sort_values('mean_test_score', ascending = False)
# Mejores hiperparámetros
# ==============================================================================
print("-----------------------------------")
print("Mejores hiperparámetros encontrados")
print("-----------------------------------")
print(grid.best_params_, ":", grid.best_score_, grid.scoring)
# Mejor modelo encontrado
# ==============================================================================
modelo = grid.best_estimator_
Predicción intervalos¶
Una vez conocidos los parámetros de la distribución condicional, se pueden calcular los cuantiles del índice de masa corporal esperados para cada edad.
# 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(datos.age), stop=max(datos.age), num=2000)
predicciones_dist = modelo.pred_dist(X=grid_X.reshape(-1, 1))
pred_cuantil_10 = [stats.lognorm.ppf(q=0.1, **dist.params) for dist in predicciones_dist]
pred_cuantil_90 = [stats.lognorm.ppf(q=0.9, **dist.params) for dist in predicciones_dist]
df_intervalos = pd.DataFrame({
'age': grid_X,
'pred_cuantil_10': pred_cuantil_10,
'pred_cuantil_90': pred_cuantil_90
})
# Gráfico intervalo cuantílico 10%-90%
# ==============================================================================
fig, axs = plt.subplots(ncols=2, nrows=1, figsize=(10, 4.5), sharey=True)
datos[datos.age <= 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[0]
)
df_intervalos[df_intervalos.age <= 2.5].plot(
x = 'age',
y = 'pred_cuantil_10',
c = 'red',
kind = "line",
ax = axs[0]
)
df_intervalos[df_intervalos.age <= 2.5].plot(
x = 'age',
y = 'pred_cuantil_90',
c = 'red',
kind = "line",
ax = axs[0]
)
axs[0].set_ylim([10, 30])
axs[0].set_title('Edad <= 2.5 años')
datos[datos.age > 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[1]
)
df_intervalos[df_intervalos.age > 2.5].plot(
x = 'age',
y = 'pred_cuantil_10',
c = 'red',
kind = "line",
ax = axs[1]
)
df_intervalos[df_intervalos.age > 2.5].plot(
x = 'age',
y = 'pred_cuantil_90',
c = 'red',
kind = "line",
ax = axs[1]
)
axs[1].set_title('Edad > 2.5 años')
axs[1].get_legend().remove()
fig.tight_layout()
plt.subplots_adjust(top = 0.85)
fig.suptitle('Distribución índice masa corporal en función de la edad', fontsize = 14);
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.lognorm.ppf(q=0.1, **dist.params) for dist in predicciones_dist]
pred_cuantil_90 = [stats.lognorm.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}")
Anomalías (outliers)¶
Conocer la distribución de probabilidad resulta útil para detectar anomalías, observaciones con valores muy poco probables acorde al modelo considerado.
Dada una determianda edad, índices de masa corporal bajos son un indicativo de posibles problemas de malnutrición y, valores muy altos, son indicativos de potenciales problemas de sobrepeso. Utilizando el modelo entrenado, se identifican aquellas personas con valores de bmi muy poco probables para la edad que tienen.
# Predicción intervalos del 98% para las observaciones disponibles
# ==============================================================================
X = datos.age.to_numpy()
y = datos.bmi.to_numpy()
predicciones_dist = modelo.pred_dist(X=X.reshape(-1, 1))
pred_cuantil_01 = [stats.lognorm.ppf(q=0.01, **dist.params) for dist in predicciones_dist]
pred_cuantil_99 = [stats.lognorm.ppf(q=0.99, **dist.params) for dist in predicciones_dist]
dentro_intervalo = np.where((y >= pred_cuantil_01) & (y <= pred_cuantil_99), True, False)
df_intervalos = pd.DataFrame({
'age': X,
'bmi': y,
'pred_cuantil_01': pred_cuantil_01,
'pred_cuantil_99': pred_cuantil_99,
'dentro_intervalo': dentro_intervalo
})
df_intervalos.head()
Se resaltan en rojo aquellos individuos para los que, aocorde al modelo, tienen una probabilidad de ocurrir igual o inferior al 1%.
# Gráfico anomalías
# ==============================================================================
fig, axs = plt.subplots(ncols=2, nrows=1, figsize=(10, 4.5), sharey=True)
df_intervalos[df_intervalos.age <= 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[0]
)
df_intervalos[(df_intervalos.age <= 2.5) & (df_intervalos.dentro_intervalo == False)].plot(
x = 'age',
y = 'bmi',
c = 'red',
kind = "scatter",
label = 'anomalia',
ax = axs[0]
)
axs[0].set_ylim([10, 30])
axs[0].set_title('Edad <= 2.5 años')
df_intervalos[df_intervalos.age > 2.5].plot(
x = 'age',
y = 'bmi',
c = 'black',
kind = "scatter",
alpha = 0.1,
ax = axs[1]
)
df_intervalos[(df_intervalos.age > 2.5) & (df_intervalos.dentro_intervalo == False)].plot(
x = 'age',
y = 'bmi',
c = 'red',
kind = "scatter",
ax = axs[1]
)
axs[1].set_title('Edad > 2.5 años')
fig.tight_layout()
plt.subplots_adjust(top = 0.85)
fig.suptitle('Distribución índice masa corporal en función de la edad', fontsize = 14);
Información de sesión¶
from sinfo import sinfo
sinfo()
Bibliografía¶
NGBoost: Natural Gradient Boosting for Probabilistic Prediction by Tony Duan, Anand Avati, Daisy Yi Ding, Khanh K. Thai, Sanjay Basu, Andrew Y. Ng, Alejandro Schuler arXiv:1910.03225
NGBoost: Natural Gradient Boosting for Probabilistic Prediction
An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)
Introduction to Machine Learning with Python: A Guide for Data Scientists
Python Data Science Handbook by Jake VanderPlas
Applied Predictive Modeling by Max Kuhn and Kjell Johnson
The Elements of Statistical Learning by T.Hastie, R.Tibshirani, J.Friedman
¿Cómo citar este documento?
Gradient boosting probabilístico (NGBoost) con python por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py26-gradient-boosting-probabilistico-python.html 
¿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 contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.