Capitulo12
Sharon Rodriguez
2023-10-21
R Markdown
TSstudio El paquete TSstudio nos brinda un conjunto de herramientas de análisis descriptivo y predictivo de datos de series temporales. Eso incluye funciones de utilidad para preprocesar datos de series temporales, funciones de visualización interactiva que es complementado con el paquetes plotly y un conjunto de herramientas para evaluar los modelos de pronóstico de series temporales
library(TSstudio)data(USVSales)ts_plot(USVSales,
title = "consumo mensual de gas natural",
Ytitle = "pies cubicos por millon",
Xtitle = "Años")ts_decompose(USVSales)#Analisis estacional
USVSales_detrend <- USVSales - decompose (USVSales) $trend
ts_seasonal(USVSales_detrend, type = "box")## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observationsPodemos ver en el gráfico estacional anterior que, normalmente, el pico del año se produjo durante los meses de marzo, mayo y junio. Además, se puede ver que las ventas decaen a partir de los meses de verano y vuelven a alcanzar su punto máximo en diciembre durante las temporadas navideñas. Por otro lado, el mes de enero suele ser el mes más bajo del año en términos de ventas.
##Analisis de Correlación
ts_lags(USVSales, lags = c(12, 24,36))La relación de la serie con el primer y segundo rezago estacional tiene una fuerte relación lineal, como se muestra en el gráfico de rezagos anterior.
##Analisis exploratorio
df <- ts_to_prophet (window(USVSales, start = c(2010,1)))
names (df) <- c ("date", "y")
head (df)## date y
## 1 2010-01-01 712.469
## 2 2010-02-01 793.362
## 3 2010-03-01 1083.953
## 4 2010-04-01 997.334
## 5 2010-05-01 1117.570
## 6 2010-06-01 1000.455ts_plot(df,
title = "total de vehiculos vendidos en US (subset)",
Ytitle = "uninades por mil",
Xtitle = "años")#Ingenieria de caracteristicas.
La ingeniería de características juega un papel fundamental cuando se modela con algoritmos de aprendizaje automático. Nuestro siguiente paso, basado en las observaciones anteriores, es crear nuevas características que se puedan utilizar como entrada informativa para el modelo. En el contexto de la previsión de series temporales, aquí hay algunos ejemplos de posibles nuevas características que se pueden crear a partir de la propia serie.
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, uniondf <- df %>% mutate(month = factor (month(date, label = TRUE), ordered = FALSE),
lag12 = lag (y, n = 12 )) %>%
filter(!is.na(lag12))df$trend <- 1: nrow(df)
df$trend_sqr <- df$trend ^ 2
str (df)## 'data.frame': 108 obs. of 6 variables:
## $ date : Date, format: "2011-01-01" "2011-02-01" ...
## $ y : num 836 1007 1277 1174 1081 ...
## $ month : Factor w/ 12 levels "ene","feb","mar",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ lag12 : num 712 793 1084 997 1118 ...
## $ trend : int 1 2 3 4 5 6 7 8 9 10 ...
## $ trend_sqr: num 1 4 9 16 25 36 49 64 81 100 ...#Entrenamiento, pruebas y evaluación de modelos.
h <- 12
train_df <- df[1:(nrow(df)-h),]
test_df <- df [(nrow(df)-h + 1) : nrow (df),]Anteriormente, la variable h representaba el horizonte de pronóstico, que, en este caso, también es igual a la longitud de la partición de prueba. Evaluaremos el rendimiento del modelo en función de la puntuación MAPE en la partición de prueba.
Además de las particiones de entrenamiento y prueba, necesitamos crear las entradas para el pronóstico en sí. Crearemos un dato. encuadre con las fechas de los siguientes 12 meses y construya el resto de características
Por último, pero no menos importante, extraeremos las últimas 12 observaciones de la serie del objeto df y las asignaremos como rezagos futuros de la serie
#Modelo de referencia.
Introducido en el Capítulo 8, Estrategias de pronóstico, el desempeño de un modelo de pronóstico debe medirse por la tasa de error, principalmente en la partición de prueba, pero también en la partición de entrenamiento. Debe evaluar el desempeño del modelo con respecto a algún modelo de referencia. En los capítulos anteriores, comparamos el pronóstico de la serie de gas estadounidense con el uso del modelo ingenuo estacional. En este capítulo, dado que utilizamos una familia de modelos de regresión de ML, tiene más sentido utilizar un modelo de regresión como punto de referencia para los modelos de ML. Por lo tanto, entrenaremos un modelo de regresión lineal de series de tiempo, que presentamos en el Capítulo 9, Pronósticos con regresión lineal. Usando las particiones de entrenamiento y prueba que creamos anteriormente, entrenemos el modelo de regresión lineal y evaluemos su rendimiento con las particiones de prueba.
lr <- lm(y ~ month + lag12 + trend + trend_sqr, data = train_df)Usaremos la función summary para revisar los detalles del modelo.
summary(lr)##
## Call:
## lm(formula = y ~ month + lag12 + trend + trend_sqr, data = train_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -146.625 -38.997 0.111 39.196 112.577
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 542.93505 72.59490 7.479 7.91e-11 ***
## monthfeb 112.73160 34.16141 3.300 0.001439 **
## monthmar 299.20932 54.24042 5.516 4.03e-07 ***
## monthabr 182.52406 42.53129 4.292 4.88e-05 ***
## monthmay 268.75603 51.28464 5.240 1.24e-06 ***
## monthjun 224.66897 44.26374 5.076 2.41e-06 ***
## monthjul 177.88564 42.21898 4.213 6.49e-05 ***
## monthago 241.63260 47.00693 5.140 1.86e-06 ***
## monthsept 152.99058 37.04199 4.130 8.76e-05 ***
## monthoct 125.16484 35.04896 3.571 0.000601 ***
## monthnov 127.97288 34.18772 3.743 0.000338 ***
## monthdic 278.67994 51.09552 5.454 5.21e-07 ***
## lag12 0.33906 0.10738 3.158 0.002236 **
## trend 7.73667 1.72415 4.487 2.36e-05 ***
## trend_sqr -0.05587 0.01221 -4.576 1.69e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 59.6 on 81 degrees of freedom
## Multiple R-squared: 0.9198, Adjusted R-squared: 0.9059
## F-statistic: 66.36 on 14 and 81 DF, p-value: < 2.2e-16A continuación, predeciremos los valores correspondientes de la serie en la partición de prueba con la función predict usando test_df como entrada.
test_df$yhat <- predict(lr, newdata = test_df)ahora podemos evaluar el rendimiento del modelo en la partición de prueba.
mape_lr <- mean(abs(test_df$y - test_df$yhat) / test_df$y)
mape_lr## [1] 0.03594578Por lo tanto,la puntuación MAPE del modelo de pronostico de regresion es de 3.0% usaremos esto para comparar el rendimiento de los modelos ML.
#Iniciar un grupo de h2o.
El paquete h2o se basa en el uso de computación distribuida y paralela para acelerar el tiempo de computo y poder escalar para big data. todo esto se realiza en clusteres de procesamiento distribuido en memoria (basado en la RAM interna de la computadora) o en paralelo (por ejemplo, AWS, google cloud, etc) Por lo tanto, cargaremos el paquete y luego configuraremos el cluster en memoria con h2o.
library(h2o)##
## ----------------------------------------------------------------------
##
## Your next step is to start H2O:
## > h2o.init()
##
## For H2O package documentation, ask for help:
## > ??h2o
##
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit https://docs.h2o.ai
##
## ----------------------------------------------------------------------##
## Attaching package: 'h2o'## The following objects are masked from 'package:lubridate':
##
## day, hour, month, week, year## The following objects are masked from 'package:stats':
##
## cor, sd, var## The following objects are masked from 'package:base':
##
## %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
## colnames<-, ifelse, is.character, is.factor, is.numeric, log,
## log10, log1p, log2, round, signif, trunch2o.init(max_mem_size = "16G")## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 4 minutes 43 seconds
## H2O cluster timezone: America/Bogota
## H2O data parsing timezone: UTC
## H2O cluster version: 3.42.0.2
## H2O cluster version age: 3 months and 4 days
## H2O cluster name: H2O_started_from_R_57301_asu791
## H2O cluster total nodes: 1
## H2O cluster total memory: 13.90 GB
## H2O cluster total cores: 8
## H2O cluster allowed cores: 8
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## R Version: R version 4.3.1 (2023-06-16 ucrt)h2o.init le permite configurar el tamaño de la memoria del cluster con el argumento max_mem_size. El resultado de la función, como se muestra en el codigo anterior, proporciona información sobre la configuración del cluster.
Cualquier dato que se utilice durante el proceso de entranamiento y prueba de los modelos mediante el paquete h2o debe cargarse en el propio cluster. La función as.h2o nos permite transformar cualquier objeto data.frame en un cluster h2o.
train_h <- as.h2o(train_df)##
|
| | 0%
|
|======================================================================| 100%test_h <- as.h2o(test_df)##
|
| | 0%
|
|======================================================================| 100%Ademas, tranformaremos el objeto predict_df (los valores futuros de las entradas de la serie) en un objeto h2o, que se utilizará para generar, mas adelante en este capitulo, el pronostico final.
Etiquetaremos los nombres de las variables dependientes e independientes.
x <- c("month","lag12", "trend", "trend_sqr")
y <- "y"#Entrenando un modelo de ML
El paquete h2o proporciona un conjunto de herramientas para entrenar y probar modelos de ML. los enfoques de entranamiento de modelos mas comunes son los siguientes. Alo largo de este capitulo, utiliozaremos el enfoque de validación cruzada (CV) para entrenar estos modelos.
#Previsión con el modelo random forest.
Ahora que hemos preparado los datos, creado nuevas funciones y lanzado un clúster de agua, estamos listos para construir nuestro primer modelo de pronóstico con el algoritmo Random Forest (RF). El algoritmo RF es uno de los modelos ML más populares y se puede utilizar tanto para problemas de clasificación como de regresión. En pocas palabras, el algoritmo de RF se basa en un conjunto de múltiples modelos de árboles.
Comenzaremos con un modelo de RF simplista utilizando 500 árboles y entrenamiento CV de 5 carpetas. Además, agregaremos un criterio de parada para evitar que el modelo se ajuste al modelo mientras no haya cambios significativos en el rendimiento del modelo. En este caso, estableceremos la métrica de parada como RMSE, la tolerancia de parada como 0,0001 y los redondeos de parada a 10.
El siguiente grafico demuestra el proceso de aprendizaje del modelo en función del del números de arboles.
rf_md <- h2o.randomForest(training_frame = train_h,
nfolds = 5,
x = x,
y = y,
ntrees = 500,
stopping_rounds = 10,
stopping_metric = "RMSE",
score_each_iteration = TRUE,
stopping_tolerance = 0.0001,
seed = 1234)##
|
| | 0%
|
|======================================================================| 100%La función h2o.randomForest devuelve un objeto con información sobre la configuración de losparametros del modelo y su rendimiento en el conjunto de entrenamiento y validación, si se usa. Podemos ver la contribución de las entradas del modelo con el agua. Fución varim_plot. Esta función devuelve un grafico con la clasificación de la contribución de las variables de entrada al rendimiento del modelo utilizando una escala entre 0 y 1, como se muestra en el siguiente codigo:
h2o.varimp_plot(rf_md)
Como podemos ver en el grafico de importancia de la variable anterior, la variable de retraso, lag12, es la mas importate para el modelo. Esto no deberia ser una sorpresa ya que vimos la fuerte relación entre la serie y su defase estacional en el analisis de correlación. Inmediatamente despúes de esto, las variables mas importantes son trend_sqr, mes y tendencia.
La salis del modelo contiene ademas del modelo en si, información sobre el rendimiento y los parametros del modelo. repasemos el resumen del modelo:
rf_md@model$model_summary## Model Summary:
## number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1 41 41 31587 8
## max_depth mean_depth min_leaves max_leaves mean_leaves
## 1 12 10.04878 45 66 56.70732Podemos ver que utilizamos solo 35 arboles de los 500 establecidos por el argumento ntrees. Esto se debe a los parámetros de parada que se utilizaron en el modelo. El siguiente grafico demuestra el proceso de aprendizaje del modelo en función del número de arboles:
library(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layouttree_score <- rf_md@model$scoring_history$training_rmse
plot_ly(x = seq_along(tree_score), y = tree_score,
type = "scatter", mode = "line") %>%
layout(title = "The trained model score history",
yaxis = list(title = "RMSE"),
xaxis = list(tltle = "Num. of trees"))## Warning: Ignoring 1 observationsPor ultimo pero no meos importante, medimos el rendimiento del modelo en la participación de prueba. usaremos la función h2o.predict para predecir los valores correspondientes de la serie en la partición de prueba.
test_h$pred_rf <- h2o.predict(rf_md, test_h)##
|
| | 0%
|
|======================================================================| 100%A continuación, podemos transferir el marco de datos de agua a un archivo, obketo de marco con la función as.data.frame
test_1 <- as.data.frame(test_h)Ahora podemos calcular la puntuación MAPE del modelo Rf en la partición de prueba:
mape_rf <- mean(abs(test_1$y - test_1$pred_rf) / test_1$y)
mape_rf## [1] 0.04647164Como se puede ver en la puntuación de error el modelo RF con su configuración predeterminado tuvo una tasa de error mas alta que nuestro modelo de referencia, es decir, el modelo de regresión lineal con una puntuación Mape de de 4.9% en comparación con el 3.0%.
Demostraremos un enfoque de búsqueda de cuadrícula para ajustar un modelo de ML con la función h2o.grid. Más adelante en este capítulo, veremos un enfoque basado en algoritmos para ajustar un modelo N ML con la función h2o.automl.
La función h2o.grid le permite establecer un conjunto de valores para algunos parámetros seleccionados y probar su rendimiento en el modelo para identificar los parámetros de ajuste que optimizan el rendimiento del modelo.
hyper_params_rf <- list(mtries = c(2, 3, 4),
sample_rate = c(0.632, 0.8, 0.95),
col_sample_rate_per_tree = c(0.5, 0.9, 1.0),
max_depth = c(seq(1, 30, 3)), min_rows = c(1, 2, 5, 10))cuantos mas parametros agreguemos o definamos para una amplia gama de valores, mayor será la combinación de busqueda posible.
Tendria sentido utilizar el metodo cartesiano cuando el número de combinaciones de busqueda sea bastante pequeño o cuando no este limitado en el tiempo de busqueda. Por otro lado, para una gran cantidad de combinaciones de busqueda o cuando tenemos limitación de tiempo se recomienda utilizar la busqueda alaeatoria. Por razones de ficiencia, establecemos una busqueda aleatoria y restringiremos el tiempo de busqueda a 20 minutos con max_rutime_sec. Usaremos la misma metrica de parada que usamos anteriormente.
search_criteria_rf <- list(strategy = "randomDiscrete",
stopping_metric = "rmse",
stopping_tolerance = 0.0001,
stopping_rounds = 10,
max_runtime_secs = 60 * 20) #Pronosticar con el modelo GBM
El algorimo GBM es otro modelo basado en conjuntos y arboles. utiliza el enfoque de impulso para entrenar diferentes subconjuntos de datos y repite el entrenamiento de subconjunto que tenia el modelo con una alta tasa de error. Esto permite que el modelo aprenda de los errores pasados y mejore el poder predictivo del modelo.
El siguiente ejemplo demuestra el uso de la función h2o.gbm para entrenar el model CBM con los mismos datos de entrada que usamos anteriormente.
De manera similar el modelo RF, podemos revisar el rango de importancia de las variables del modelo con la función h2o.varimp_plot:
para RF, el modelo GBM clasifica la variable lag12 como la mas importante para el modelo.Probaremos el rendimiento del modelo n el conjunto de prueba.
El modelo GBM obtuvo una tasa de error del 3.9%, en comparación con los otros modelos que hemos probado que uno fue del 3.0% y otro del 4.9%
Visualicemos los resultados y comparemos la predicción con la lunea de base y la real predicción:
library(plotly)
plot_ly(data = test_1) %>%
add_lines(x = ~ date, y = ~y, name = "Actual") %>%
add_lines(x = date, y = ~ yhat, name = "Linear Regression", line = list(dash = "dot")) %>%
add_lines(x = ~ date, y= ~ pred_rf, name = "random Forest", line = list(dash = "dash")) %>%
add_lines(x = ~ date, y = ~ grid, name = "random Forest (grid)", line = list(dash = "dash")) %>%
layout(tittle = "Total Vehiculos Vendidos - Actual vs. prediction (random Forest)",
yaxis = list(title = "Thousands of units"),
xaxis = list(title = "Month"))## Warning: 'layout' objects don't have these attributes: 'tittle'
## Valid attributes include:
## '_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'smith', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'#Previsión con el modelo AutoMl
Hasta ahora, en este capítulo, hemos analizado dos enfoques de modelado: el primero utiliza la configuración predeterminada del algoritmo con los modelos RF y GBM, y el segundo aplica una búsqueda de cuadrícula con el modelo RF. Echemos un vistazo a un tercer enfoque para ajustar los modelos de ML. Aquí usaremos la función h2o.automl. El agua. La función automl proporciona un enfoque automatizado para entrenar, ajustar y probar múltiples algoritmos de ML antes de seleccionar el modelo que funcionó mejor según la evaluación del modelo. Utiliza algoritmos como RF, GBM, DL y otros que utilizan diferentes enfoques de ajuste.
De manera similar, la función h2o.grid puede aplicar cualquiera de los enfoques de entrenamiento (CV, entrenamiento con validación, etc.) durante el proceso de entrenamiento de los modelos. Usemos la misma entrada que antes y entrenemos el modelo de pronóstico.
El modelo líder de producción de agua. autom1 logró una puntuación MAPE del 4.2%. Aunque los modelos anteriores lograron una puntuación más baja todavía proporciona un aumento significativo con respecto al modelo base. Visualicemos la predicción del modelo en las particiones de prueba con respecto a la predicción real y de referencia
`
La principal ventaja de la función h2o.automl es que puede ampliarse y tener varias series para pronosticar con una mínima intervención por parte del usuario. Esto, por supuesto, conlleva potencia de cálculo adicional y coste de tiempo.
#Seleccionando el modelo final.
Dado que todos estos modelos lograron mejores resultados que el modelo de referencia, podemos descartar el modelo de referencia. Además, dado que el segundo modelo de RF (con búsqueda de cuadrícula) logró mejores resultados que el primero, no tiene sentido conservar el primer modelo. Esto nos deja con tres modelos de pronóstico, es decir, RF (con búsqueda de cuadrícula), GBM y AutoML. Generalmente, dado que el modelo GBM logró los mejores resultados MAPE, lo seleccionaremos. Sin embargo, siempre es bueno trazar el pronóstico real y comprobar cómo se ve el pronóstico real. Antes de trazar los resultados, usemos estos tres modelos para pronosticar los próximos 12 meses usando los datos. pronóstico de marco que creamos en las secciones Pronóstico con el modelo Random Forest y Pronóstico con el modelo GBM.