XGBoost para predicción en Series Temporales

Contenido¶

• Modelado
• Estructuración serie temporal como problema de regresión
• Conjunto de datos
  • Descripción variables
  • Preparación datos
  • Distribución datos
• Modelo XGBoost
  • Datos Train y Test
  • Entrenamiento modelo
    • Hiperparámetros
    • Evaluación desempeño modelo
    • Importancia de los predictores

Modelado¶

Existen multitud de métodos de previsión de la demanda para cualquier sector de bienes y servicios.Para nuestro propósito realizamos diferentes implementaciones sobre el mismo objetivo en el pronóstico.Por una lado, el modelo ARIMA (AutoRegresive Integrates Moving Average) o modelo autorregresivo integrado de media móvil en español, y por otro lado en este caso, la implementación del método de ensemble por medio de la librería XGboost.

Aprovecharé para resaltar algunas características por las que se hace uso extendido de la librería.

El primer modelo mencionado es conocido por su uso tradicional y cuyos resultados han sido efectivos. El segundo modelo incorpora variables adicionales además del componente temporal, creando un ajuste multivariante. En el caso de ARIMA se trata de un modelo univariante, pues no toma en consideración otras variables en el pronóstico.

Para mayor riqueza y sobretodo precisión, en un entorno de producción real, una de las etapas o ciclos previos al modelaje de datos,sería la validación de datos y continuidad en el registro temporal, o por ejemplo cálculos numéricos del error (se realizan cálculos sobre las ventas y se producen pérdidas por redondeo) entre otros aspectos.

Estructuración serie temporal como problema de regresión¶

En el caso de las series temporales,a comparación del planteamiento de los datos (de por ejemplo clasificación) como puntos independientes en los modelos de aprendizaje supervisado, hay que tratar la componente temporal. La dimensión temporal añade un orden explícito de los datos que ha de preservarse porque provee información adicional/importante a los algoritmos de aprendizaje.

Para reflejar la tendencia de ventas en el tiempo de los datos, por un lado se creará una variable con la consecución de meses como números. Por ejemplo, la fecha 2012-01-03 tendrá el valor 0 como inicial y 2020-01-15 será el mes 96.

Por otro lado, se obtiene la media móvil simple ya que ofrece una visión suavizada de la serie.Promediando varios valores se elimina parte de los movimientos irregulares de la serie.

Los modelos de previsión de ventas en aplicaciones industriales han ido desarrollándose entorno a su componente temporal,por el registro histórico que se tiene.En la práctica actualmente, y por comparativa, los modelos modernos obtienen mejor resultados. Se entiende que la predicción de ventas es más bien un problema de regresión, y no deteminante únicamente por temporalidad.

Conjunto de datos¶

El conjunto de datos (+3 Millones de instancias) contiene las compras al por mayor de licores en el Estado de Lowa por minoristas para la venta por unidades desde enero del año 2012 hasta el 2020. Ya que estamos trabajando sobre la memoria local y el conjunto de datos es relativamente grande, implementaremos una solución multiproceso que reduce el tiempo de ejecución respecto a la solución secuencial.

El Estado de Lowa controla la distribución al por mayor de licores destinados a la venta al por menor, lo que significa que el conjunto de datos ofrece una vista completa del sector de la venta de licores en el Estado.

En nuestro caso, el enfoque será a partir de los datos de un proveedor seleccionado, para el cual los datos sí son continuos en el tiempo y el tamaño del conjunto es considerable.

Descripción variables¶

Las variables de nuestro interés son las siguientes:

Variable	Descripción
Invoice/Item Number	Concatenación de número de factura y línea asociada con el pedido del licor
Date	Fecha de pedido
County Number	Número del condado donde la tienda que ha realizado el pedido se localiza
Category Name	Nombre de la categoría del licor pedido
Bottle Volume (ml)	Volumen de cada botella de licor pedida en milimetros
State Bottle Cost	Total que recibe la División de Bebidas Alcohólicas (administración) por cada botella de licor pedida
State Bottle Retail	Total que paga cada tienda por cada botella de licor pedida
Bottles Sold	Número de botellas de licor pedidas por la tienda
Sale (Dollars)	Coste total de pedido de licores
Volume Sold (Liters)	Volumen total de licor pedido en litros

Preparación datos¶

En la implementación multiproceso se han realizado las siguientes operaciones,aunque en esta publicación se realizan adecuaciones adicionales que se detallan más abajo de las variables para el modelo.

 * Renombrar columnas
 * Encontrar datos duplicados
 * Revisión valores perdidos y ceros
 * Tratamiento de tipos de datos y valores únicos
   * Agrupacion de datos en grupos reducidos
 * Exploración distribución variables

A modo de síntesís, respecto al tratamiento de datos realizado previamente:

• Se realiza conversión de la fecha
• Sobre la categoría de licores, se ha detectado que en general que los items
  pertencen a variedad de whiskeys, vodkas o rones. Mediante búsqueda por cadenas de
  texto se realiza una nueva clasificación.
• En cuanto al tamaño de las botellas, hay más de 30 tamaños diferentes (volume by ml). Sin embargo, la mayoría de botellas vendidas son de 750ml, 1000ml o 17500ml.
  Para manejar mejor la variables, se revierte la columna indicando si el tamaño de la botella de licor es de 750ml, 1000ml, 1750ml o de otro tamaño.

• El porcentaje de valores peridos es insignificante (1%).

A fin de ajustar los datos al modelo actual:

• Se descarta la variable invoice_line referente al número de factura tras comprobar la duplicidad de los registros

• Para los valores perdidos de county_no referente al número para identificar el condado, hemos asignado valor 'UNKNOWN' ya que el tratamiento será categórica y de esta forma no perdemos número de registros

• Se realiza cálculo del número consecutivo de los meses y media móvil mensual(indexado) por categoría y tipo de tamaño de botellas y se vuelve a reindexar el dataframe tras elimnar el campo Date.

In [ ]:

# Importar librerias
import pandas as pd
import numpy as np
from dateutil.relativedelta import relativedelta # intervalos de tiempo
import sklearn
from sklearn.metrics import explained_variance_score,  mean_squared_error, r2_score 
import xgboost as xgb
from xgboost import XGBRegressor

import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px

import warnings
warnings.filterwarnings('ignore')

In [ ]:

class vendor:
  def __init__(self,dir_name):
    self.dir_name = dir_name

  
  def read_file(self):
    df = pd.read_csv(self.dir_name,
                     parse_dates=['Date'],
                     dtype={'county_no': 'object'}).sort_values(by='Date')
    
    # Revision valores duplicados
    duplicate = df[df.duplicated()]
    print("Exiten {} instancias duplicadas \n\n".format(len(duplicate)))
    print("Tamaño del conjunto de datos original es de {} observaciones y {} variables \n\n".format(df.shape[0],
                                                                                  df.shape[1]))
    print("Columnas de conjunto de variables:\n{}\n\n".format(list(df.columns)))
    
    return df

  
  def extract_innecesary_features(self):
    df = vendor.read_file(self)
    # Eliminar invoice_line
    df = df.drop(['invoice_line_no'], axis=1)
    
    # el tratamiento del numero de condado es categorico. asignamos valor unknown en los registros con valor nan
    df['county_no'].fillna('UNKNOWN', inplace=True)

    # Columnas valor incial 0 para numero consecutivo de meses y media movil simple
    df['date_block_num'] = 0
    df['SMA'] = 0

    # Agrupacion mensual por nombre de categoria y volumen de botella para crear nuevas variables
    df = df.set_index('Date').groupby(['category_name','bottle_volume_ml']).resample('M').mean()

    print("\n\t Datos tras la agrupcion mensual por categoria y tipo de botella\n\t")
    display(df.head(5))
    
    return df
   

  def feature_engineering(self):
    df = vendor.extract_innecesary_features(self)

    # fecha inicial para calcular intervalos de tiempo en meses
    first = []

    # Para cada categoria, tipo de tamaño de botella y sus correspondientes meses,
    # se realiza el calculo de meses consecutivos siendo el valor del  primer mes 0
    for cat in list(df.index.get_level_values('category_name').unique()):
      for bottle in list(df.index.get_level_values('bottle_volume_ml').unique()):
        for a,date in df.xs([cat,bottle]).groupby(level=0):
          date1 = first.append(pd.to_datetime(date.index.strftime('%Y-%m-%d %H:%M:%S')[0]))
          date2 = pd.to_datetime(date.index.strftime('%Y-%m-%d %H:%M:%S')[0])
          r = relativedelta(date2, first[0])
          df.xs([cat, bottle, \
                 pd.to_datetime(date.index.strftime('%Y-%m-%d %H:%M:%S')[0]) ]).date_block_num = int(r.years*12 + r.months)

    # calculo de media movil para cada tipo de categoria y tipo de tamaño de botella
    for cat in list(df.index.get_level_values('category_name').unique()):
      for bottle in list(df.index.get_level_values('bottle_volume_ml').unique()):
        df.xs([cat, bottle])['SMA'] = df.xs([cat,bottle]).sale_dollars.rolling(12).mean()
    

    return df
      



  def missing_values(self):
    df = vendor.feature_engineering(self)
    
    # Reindexar dataframe
    df = df.reset_index().drop(['Date'], axis=1)
    
    # Nan por valores constantes 0
    df.SMA.fillna(0, inplace=True)

    print("\n\t Datos tras crear columna con numero consecutivo de meses y media movil con ventada=12 para las categorias,\n\
    reindexar \n\t")
    display(df.head(12))
    return df


  def one_hot_encoder(self):
    df = vendor.missing_values(self)
    # cambiar nombre variable target a 'y'
    df = df.rename({'sale_dollars': 'y'}, axis=1)
    #dummies variables
    dfDummies = pd.get_dummies(df[['category_name', 'bottle_volume_ml']], prefix=['category', 'bottle'])
    vendor_df = pd.concat([df,dfDummies], axis=1).drop(['category_name', 'bottle_volume_ml'],axis=1)
    print("\n\n Datos preparados para modelar \n\n")
    display(vendor_df.head())
    return vendor_df
    

if __name__ == '__main__':
  dir_name = '/content/drive/My Drive/prepared_data_liquor_sales/data_arima_cleaned.csv'
  vendor_data = vendor(dir_name)
  vendor_data.one_hot_encoder()

Exiten 0 instancias duplicadas 


Tamaño del conjunto de datos original es de 3054461 observaciones y 10 variables 


Columnas de conjunto de variables:
['Date', 'invoice_line_no', 'county_no', 'category_name', 'bottle_volume_ml', 'state_bottle_cost', 'state_bottle_retail', 'sale_bottles', 'sale_dollars', 'sale_liters']



	 Datos tras la agrupcion mensual por categoria y tipo de botella
	

			state_bottle_cost	state_bottle_retail	sale_bottles	sale_dollars	sale_liters	date_block_num	SMA
category_name	bottle_volume_ml	Date
GIN	1000ml	2012-01-31	14.200233	21.294884	9.472868	198.284264	9.472868	0.0	0.0
		2012-02-29	14.115615	21.167907	10.083056	213.292226	10.083056	0.0	0.0
		2012-03-31	14.309792	21.459288	9.596439	205.823264	9.596439	0.0	0.0
		2012-04-30	14.181429	21.266667	13.226190	282.651012	13.226190	0.0	0.0
		2012-05-31	14.111724	21.162069	10.184729	211.850985	10.184729	0.0	0.0

	 Datos tras crear columna con numero consecutivo de meses y media movil con ventada=12 para las categorias,
    reindexar 
	

	category_name	bottle_volume_ml	state_bottle_cost	state_bottle_retail	sale_bottles	sale_dollars	sale_liters	date_block_num	SMA
0	GIN	1000ml	14.200233	21.294884	9.472868	198.284264	9.472868	0.0	0.000000
1	GIN	1000ml	14.115615	21.167907	10.083056	213.292226	10.083056	1.0	0.000000
2	GIN	1000ml	14.309792	21.459288	9.596439	205.823264	9.596439	2.0	0.000000
3	GIN	1000ml	14.181429	21.266667	13.226190	282.651012	13.226190	3.0	0.000000
4	GIN	1000ml	14.111724	21.162069	10.184729	211.850985	10.184729	4.0	0.000000
5	GIN	1000ml	14.522047	21.777795	16.498688	362.163097	16.498688	5.0	0.000000
6	GIN	1000ml	14.154918	21.226885	10.348946	214.891850	10.348946	6.0	0.000000
7	GIN	1000ml	14.078031	21.111509	9.946292	206.942660	9.946292	7.0	0.000000
8	GIN	1000ml	13.970000	20.949399	9.908127	203.917314	9.908127	8.0	0.000000
9	GIN	1000ml	14.148017	21.216529	10.030303	211.653223	10.030303	9.0	0.000000
10	GIN	1000ml	14.365318	21.542609	15.578595	333.035518	15.578595	10.0	0.000000
11	GIN	1000ml	14.228668	21.337554	9.535326	202.618098	9.535326	11.0	237.260293

 Datos preparados para modelar 

	state_bottle_cost	state_bottle_retail	sale_bottles	y	sale_liters	date_block_num	SMA	category_GIN	category_LIQUOR	category_RUM	category_UNKNOWN	category_VODKA	category_WHISKEY	bottle_1000ml	bottle_1750ml	bottle_750ml	bottle_other_ml
0	14.200233	21.294884	9.472868	198.284264	9.472868	0.0	0.0	1	0	0	0	0	0	1	0	0	0
1	14.115615	21.167907	10.083056	213.292226	10.083056	1.0	0.0	1	0	0	0	0	0	1	0	0	0
2	14.309792	21.459288	9.596439	205.823264	9.596439	2.0	0.0	1	0	0	0	0	0	1	0	0	0
3	14.181429	21.266667	13.226190	282.651012	13.226190	3.0	0.0	1	0	0	0	0	0	1	0	0	0
4	14.111724	21.162069	10.184729	211.850985	10.184729	4.0	0.0	1	0	0	0	0	0	1	0	0	0

Distribución datos¶

Variables categóricas

En un análisis previo en la preparación de datos se han obtenido los siguientes detalles:

• Predominan la venta de licores en la categoría de whiskey y vodka, que junto al ron conforman alrededor del 80% de las ventas. En cuanto al volumen de las botellas, cerca del 50% son de tamaño 750ml.

Variables numércias

De la misma manera, en cuanto a las variables numéricas:

• Lejos de tener una distribución uniforme, lo que vemos son dispersidad de valores que se podrían considerar atípicos, ya que en todos los casos, la mayor frecuencia se concentra un un determinado rango. Descartamos optar por segmentos de igual rango por la distribución.

In [ ]:

variables = ['state_bottle_cost', 	'state_bottle_retail' 	,'sale_bottles' ,	'y' ,	'sale_liters']
g = sns.PairGrid(df, vars=variables, diag_sharey=False, corner=True)
g.map_lower(sns.scatterplot)
g.map_diag(sns.kdeplot)
g.add_legend()

Out[ ]:

<seaborn.axisgrid.PairGrid at 0x7f7997fcd750>

Modelo XGBoost¶

Recalcando lo mencionado en la introducción, XGBoost es un framework popular para la implementación de algoritmos de boosting. Lo es, entre otros aspectos, por organizar y ejecutar el trabajo en computación paralela o la optimización en la gestión del almacenaje de datos.

Boosting es un método de ensemble que combina múltiples modelos predictivos para lograr un equilibrio entre bias (cómo de bien se aproxima el modelo a la relación real entre las variable) y varianza (cuánto varía el modelo dependiendo de la muestra empleada en el entrenamiento).

Consiste en ajustar secuencialmente múltiples modelos sencillos, llamados weak learners, de forma que cada modelo aprende de los errores del anterior.

$\blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft$ BREAK $\blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright$

En la implementación del modelo, me ha surgido la duda en qué se diferencia este método de ensemble y el funcionamiento de una red neuronal en el deep learning.

En ambos casos se registra el descenso conseguido a cada división (capa en el caso de la red neuronal), pero en el modelado de los datos, en la técnica boosting es posible identificar la influencia que tiene cada predictor sobre el modelo. En comparación con las redes neuronales, el resultado del peso del aprendizaje y bias se computa para todo el sistema.

$\blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft \blacktriangleleft$ CONTINUE $\blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright \blacktriangleright$

Problema regresión¶

La idea es construir varios árboles $T^K$. Cada árbol producirá un salida numérica para cada variable-vector $x$ que denotaremos $f^{k} (x)$.

En el ensemble de árboles cada nuevo modelo empleará la información del modelo anterior para aprender de sus errores, mejorando el resultado de la predicción $y(x)$ de la variable dependiente iteración a iteración.

Por lo que tendremos una secuencia de predicciones:

\begin{align} y^{(0)}(x) = \frac{1}{2} \\ y^{(1)}(x) = y^{(0)}(x) + \epsilon f^{(1)}(x) \\ y^{(k)}(x) = y^{(k-1)}(x) + \epsilon f^{(k)}(x)= \frac{1}{2} + \epsilon \sum_{a=1}{k}f^{a}(x) \end{align}

Datos Train y Test¶

Para dividir los datos de entrenamiento no recogemos los datos transversales porque se trata de una serie temporal dependiente, lo datos anteriores que se consideran predecirán los datos futuros.

Dado que tratamos de construir un modelo de predicción mensual, reservaremos los datos del último mes (96), entrenándolo con los datos de los meses previos

In [ ]:

X_train, X_test = vendor[vendor['date_block_num']<96], vendor[vendor['date_block_num']==96]
X_train['y'], X_test['y'] = X_train['y'], X_test['y']
    
y_train, y_test = X_train['y'], X_test['y']
del X_train['y']
del X_test['y']

Entrenamiento modelo¶

Los algoritmos de boosting se caracterizan por tener una cantidad considerable de hiperparámetros, cuyo valor óptimo se identifica mediante validación cruzada. En el argumento eval_set se proveen los datos de entrenamiento y test para la validación cruzada.

Hiperparámetros¶

• Número de weack learners o número de iteraciones La técnica boosting puede sufrir overfitting si el valor es excesivamente alto. Se pasa en el argumento n_estimators

• Learning rate $(\eta)$. Controla el ritmo al que aprende el modelo. Suelen recomendarse valores de 0.01 o 0.001. Usaremos el primer valor, lo que significa que más árboles se necesitan para alcanzar buenos resultados pero menor es el riesgo de overfitting. El parámetro $\eta$ que se especifique pasa por el término que regula la función de coste en el problema de regresión. Se pasa en el argumento learning_rate

• El número de divisiones viene dado a partir de $\gamma$. Cuanto mayor sea el valor de gamma, más conservador será el algoritmo (menor margen de error ya que los nodos nos variarán) Se pasa en el argumento min_split_loss

• Lasso. Método de regularización en la selección de variables.Penaliza la suma del valor absolutos de los coeficientes de regresión. A esta penalización se le conoce como l1 y tiene el efecto de forzar a que los coeficientes de los predictores tiendan a cero. Dado que un predictor con coeficiente de regresión cero no influye en el modelo, lasso consigue excluir los predictores menos relevantes.

No se incluye el desarrollo matemático completo del modelaje por extensión.

In [ ]:

class XGBoost:
  def __init__(self, vendor_df):
    self.vendor_df = vendor_df


  def train_test_split(self):
    """ Separar datos entrenamiento y test """
    vendor = vendor_df.dropna()
    X_train, X_test = vendor[vendor['date_block_num']<96], vendor[vendor['date_block_num']==96]
    X_train['y'], X_test['y'] = X_train['y'], X_test['y']
    y_train, y_test = X_train['y'], X_test['y']
    del X_train['y']
    del X_test['y']

    return X_train, X_test, y_train, y_test
  

  def hyper_parameters(self):
    """ Instancia objeto regresion XGBoost con los parametros indicados"""
    model =  XGBRegressor(n_estimators=100,learning_rate=0.1,reg_lambda=20,min_split_loss=0,silent=True)

    return model

  def fit_predict(self):
    """ Ajuste modelo y prediccion"""
    X_train, X_test, y_train, y_test = XGBoost.train_test_split(self)
    regr = XGBoost.hyper_parameters(self)
    
    regr.fit(X_train,
              y_train, 
              eval_metric=['rmse','mae'],
              eval_set=[(X_train, y_train), (X_test, y_test)],
              verbose=False)
    
    y_pred = regr.predict(X_test)
    feature_importance = regr.feature_importances_
    results = regr.evals_result()
    
    return y_pred, feature_importance, results

  

  def model_performance(self):
    "Evaluacion desempeño modelo"
    _, _, _, y_test = XGBoost.train_test_split(self)
    y_pred,_,_ = XGBoost.fit_predict(self)

    print('MSE on the test set data:',mean_squared_error(y_test, y_pred))
    print('R-squared on the test set data:',r2_score(y_test, y_pred))
  
  def model_metrics(self):
    """ Extrae los resultados y genera graficos """
    _,_,results = XGBoost.fit_predict(self)

    epochs = len(results['validation_0']['rmse'])
    x_axis = range(0, epochs)

    #Generar graficas
    f, ax = plt.subplots(nrows=1, ncols=2, figsize=(15,9))

    # plot RMSE
    sns.lineplot(x_axis, results['validation_0']['rmse'], label='Train', ax=ax[0])
    sns.lineplot(x_axis, results['validation_1']['rmse'], label='Test', ax=ax[0])
    ax[0].set_xlabel('Particiones')
    ax[0].set_ylabel('RMSE')
    ax[0].set_title('XGBoost RMSE')
    
    # plot MAE
    sns.lineplot(x_axis, results['validation_0']['mae'], label='Train', ax=ax[1])
    sns.lineplot(x_axis, results['validation_1']['mae'], label='Test', ax=ax[1])
    ax[1].set_xlabel('Particiones')
    ax[1].set_ylabel('MAE')
    ax[1].set_title('XGBoost MAE')
    plt.tight_layout()
    plt.show()


  def feature_impr(self):
    """ Grafico ponderacion de las variables """
    X_train, _, _, _ = XGBoost.train_test_split(self)
    _,feature_importance_reg,_= XGBoost.fit_predict(self)

    feature_importances = pd.DataFrame({'col': X_train.columns,'imp':feature_importance_reg})
    feature_importances = feature_importances.sort_values(by='imp',ascending=False)
    fig = px.bar(feature_importances,x='col',y='imp')
    fig.show()

Evaluación desempeño modelo¶

El alto valor de $R^2$ nos indica que el modelo explica casi por completo la variabilidad. Esto puede significar que probablemente estemos perdiendo información o el problema sí se ha aprendido por el modelo.

El Mean Square Test Error asociado con el árbol de regresión es de 442.65 unidades. La raíz cuadrado del Mean Square Test error es 210.39, lo que significa que las predicciones se alejan en promedio 21.04 unidades (2104 dólares) del valor real.

Por medio del gráfico de RMSE , podemos observar que los valores de Test se encuentran por debajo de los datos del entrenamiento, lo que nos indica que no hay sobre-ajuste de los datos.

En cuanto al gráfico de MAE, vemos que la distancia del promedio de la desviación de los valores obtenidos en la predicción y los valores actuales es bastante próxima entre sí

In [ ]:

print(xgbt.model_performance())
print("\n")
print(xgbt.model_metrics())

MSE on the test set data: 442.6471683436427
R-squared on the test set data: 0.9698087495653751
None

None

Importancia de los predictores¶

Se tiene el cómputo de importancia de los predictores. La variable principal refiere al total de botellas vendidas (por mes), seguido de la variable en la que hemos reunido todas aquellos tipos de volumen de botellas que en un principio eran 'inusuales' o 'poco frecuentes'.

In [2]:

print(xgbt.feature_impr())

Out[2]: