R 机器学习预测时间序列模型

2020-07-02  本文已影响0人  jamesjin63

机器学习在时间序列数据上应用

随着疫情的变化,急性传染病数据经常会随时间变化,我们通过对每天传染病的记录,就形成了时间序列数据,周期可以是天,周,月,年。目前我们经常会用到ARIMA来预测疾病在未来的变化趋势。


image.png

但是随着机器学习的广泛应用,在时间序列上,也可以采用机器学习发方法去预测,结果比传统的ARIMA EST更加快速,简洁,准确。
这次将要介绍关于的时间序列预测的Modeltime包,旨在加快模型评估,选择和预测的速度。modeltime通过将tidymodels机器学习软件包生态系统集成到简化的工作流中以进行tidyverse预测来实现此目的。modeltime结合了机器学习模型,经典模型和自动化模型等。

image.png

主要优点:
简化数据建模预测流程,包括数据建模,评估,预测及输出

  1. 预测的系统工作流程。modeltime_table(),modeltime_calibrate()modeltime_refit()
  2. 结合Tidymodels以期加入机器学习算法。如XGBoostGLMnetStanRandom Forest
  3. 改进传统时间序列模型。如arima_boost()prophet_boost()

1.数据

我们选取bike_sharing_daily时间序列数据集,其中包括自行车每日的使用数据。
这里只需要日期与当日的使用量“date” and “value”。然后可以简单绘制一下。
注意这里的时间序列是tibble格式。

install.packages("modeltime")
devtools::install_github("business-science/timetk")
library(tidyverse)
library(tidymodels)
library(modeltime)
library(timetk)   
library(lubridate)

# data
bike_transactions_tbl <- bike_sharing_daily %>%
  select(dteday, cnt) %>%
  set_names(c("date", "value")) 

bike_transactions_tbl

# plot
bike_transactions_tbl %>%
  plot_time_series(date, value, .interactive = FALSE)

2.模型构建

下面继续对数据进行处理,在建模之前,我们需要将数据分成traintest
使用time_series_split()来分割我们的数据,assess = "3 months"来确定后三个月为test数据集,cumulative = TRUE指定前面部分为train。

# Split your time series into training and testing sets
splits <- bike_transactions_tbl %>%
  time_series_split(assess = "3 months", cumulative = TRUE)
  
splits
 <Analysis/Assess/Total>
<641/90/731>
  
# plot
splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(date, value, .interactive = FALSE)

image.png
image.png

2.1经典ARIMA与Prophet

因为model time可以简化数据建模流程,包括简化参数的设置,自动模型通常是已经包含了自动化的建模方法。包括“自动ARIMA”和“自动ETS”功能以及“ Prophet”算法。主要包含三个参数设置:

结下来我们尝试建立ARIMA模型与Prophet模型

# ARIMA
model_fit_arima <- arima_reg() %>%
  set_engine("auto_arima") %>%
  fit(value ~ date, training(splits))

model_fit_arima 
parsnip model object

Fit time:  426ms 
Series: outcome 
ARIMA(0,1,3) with drift 

Coefficients:
          ma1      ma2      ma3   drift
      -0.6106  -0.1868  -0.0673  9.3169
s.e.   0.0396   0.0466   0.0398  4.6225

sigma^2 estimated as 730568:  log likelihood=-5227.22
AIC=10464.44   AICc=10464.53   BIC=10486.74

# Prophet
model_fit_prophet <- prophet_reg() %>%
  set_engine("prophet", yearly.seasonality = TRUE) %>%
  fit(value ~ date, training(splits))

model_fit_prophet
parsnip model object

Fit time:  455ms 
PROPHET Model
- growth: 'linear'
- n.changepoints: 25
- seasonality.mode: 'additive'
- extra_regressors: 0

2.2机器学习算法

机器学习模型前面设置比自动化经典模型更为复杂。通常在进行机器学习建模之前,对数据进行预处理,称之为workflow一般过程如下:

首先,我将使用配recipe()创建预处理数据的先前步骤。该过程使用“日期”列创建了我要建模的45个新的列。这些列包含了时间序列的详细信息及傅立叶变化的数据。

recipe_spec <- recipe(value ~ date, training(splits)) %>%
  step_timeseries_signature(date) %>%
  step_rm(contains("am.pm"), contains("hour"), contains("minute"),
          contains("second"), contains("xts")) %>%
  step_fourier(date, period = 365, K = 5) %>%
  step_dummy(all_nominal())

recipe_spec %>% prep() %>% juice()

有了recipe,我们可以建立机器学习模型了。
为什么需要recipe是因为在tidymodel里面,设置了建立机器学习模型的一套准则,感兴趣可以去:

机器学习模型

这里我们新建了glmnet与RF模型。
主要参数包括:

# Elastic NET model
model_spec_glmnet <- linear_reg(penalty = 0.01, mixture = 0.5) %>%
  set_engine("glmnet")

workflow_fit_glmnet <- workflow() %>%
  add_model(model_spec_glmnet) %>%
  add_recipe(recipe_spec %>% step_rm(date)) %>%
  fit(training(splits))

# Random Forest
model_spec_rf <- rand_forest(trees = 500, min_n = 50) %>%
  set_engine("randomForest")

workflow_fit_rf <- workflow() %>%
  add_model(model_spec_rf) %>%
  add_recipe(recipe_spec %>% step_rm(date)) %>%
  fit(training(splits))

2.3机器学习结合传统算法

Prophet Boost算法将Prophet与XGBoost结合使用,从而获得了两全其美的效果(即Prophet Automation + Machine Learning)

  1. 首先使用prophet对单个时间序列建模
  2. 使用通过预处理配方提供的回归数据(生成的45个新列),并使用XGBoost模型对prophet残差进行回归
model_spec_prophet_boost <- prophet_boost() %>%
  set_engine("prophet_xgboost", yearly.seasonality = TRUE) 

workflow_fit_prophet_boost <- workflow() %>%
  add_model(model_spec_prophet_boost) %>%
  add_recipe(recipe_spec) %>%
  fit(training(splits))
  
## [05:22:38] WARNING: amalgamation/../src/learner.cc:480: 
## Parameters: { validation } might not be used.
## 
##   This may not be accurate due to some parameters are only used in language bindings but
##   passed down to XGBoost core.  Or some parameters are not used but slip through this
##   verification. Please open an issue if you find above cases.

workflow_fit_prophet_boost

3.模型调整

image.png

Modeltime工作流程旨在加速模型评估和选择。
现在我们有了几个时间序列模型,让我们对其进行分析,并通过模型时间工作流程预测未来变化趋势。
Modeltime使用ID来定位我们之前建立的模型,以帮助我们识别模型。
让我们将模型添加到modeltime_table()中。
[图片上传中...(image.png-ccae46-1593703307979-0)]

model_table <- modeltime_table(
  model_fit_arima, 
  model_fit_prophet,
  workflow_fit_glmnet,
  workflow_fit_rf,
  workflow_fit_prophet_boost
) 

model_table
## # Modeltime Table
## # A tibble: 5 x 3
##   .model_id .model     .model_desc              
##       <int> <list>     <chr>                    
## 1         1 <fit[+]>   ARIMA(0,1,3) WITH DRIFT  
## 2         2 <fit[+]>   PROPHET                  
## 3         3 <workflow> GLMNET                   
## 4         4 <workflow> RANDOMFOREST             
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS

模型校准用于量化误差并估计置信区间。
我们将使用modeltime_calibrate()函数对test数据进行模型校准。
将生成两个新列(“ .type”“ .calibration_data”),其中最重要的是“ .calibration_data”
包括l测试集的实际值,拟合值和残差。

calibration_table <- model_table %>%
  modeltime_calibrate(testing(splits))

calibration_table
## # Modeltime Table
## # A tibble: 5 x 5
##   .model_id .model     .model_desc               .type .calibration_data
##       <int> <list>     <chr>                     <chr> <list>           
## 1         1 <fit[+]>   ARIMA(0,1,3) WITH DRIFT   Test  <tibble [90 × 4]>
## 2         2 <fit[+]>   PROPHET                   Test  <tibble [90 × 4]>
## 3         3 <workflow> GLMNET                    Test  <tibble [90 × 4]>
## 4         4 <workflow> RANDOMFOREST              Test  <tibble [90 × 4]>
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS Test  <tibble [90 × 4]>

4.模型预测

利用校准后的数据,我们可以对test数据进行预测。
使用modeltime_forecast()生成测试集的预测数据。
使用plot_modeltime_forecast()绘制预测结果。

calibration_table %>%
  modeltime_forecast(actual_data = bike_transactions_tbl) %>%
  plot_modeltime_forecast(.interactive = FALSE)
image.png

5.模型评估

接下来,计算预测值与实际值的差异来评估模型好坏。
使用modeltime_accuracy()生成样本外准确性指标数据。
使用table_modeltime_accuracy()生成预测结果表格

calibration_table %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(.interactive = FALSE)
image.png

后续还有对最优模型进行预测,使用全部数据作为training,预测未来1年的变化趋势。

参考

Modeltime GitHub Page - Give it a Star if you like it!
Timetk Documentation - Data wrangling, visualization, and preprocessing for time series.
Tidymodels.org - The tidymodels framework is a collection of packages for modeling and machine learning using tidyverse principles.
Introducing Modeltime

上一篇下一篇

猜你喜欢

热点阅读