Phân tích chuỗi thời gian với H2O và Timetk

Giới thiệu về H2O và Timetk trong phân tích chuỗi thời gian

h2o là một thư viện mạnh mẽ cung cấp các thuật toán học máy tiên tiến, các chỉ số hiệu suất và các hàm tiện ích, giúp việc triển khai học máy trở nên dễ dàng và hiệu quả. Một trong những lợi ích chính của h2o là khả năng triển khai trên các cụm tính toán. Từ góc độ R, h2o có bốn ứng dụng chính:

  1. Thao tác dữ liệu: Bao gồm nối, nhóm, xoay, chuyển đổi và chia tập dữ liệu thành tập huấn luyện, kiểm tra, xác thực, v.v.
  2. Thuật toán học máy: Cung cấp các thuật toán học máy giám sát và không giám sát phức tạp. Các thuật toán học máy giám sát bao gồm học sâu (mạng nơ-ron), rừng ngẫu nhiên, mô hình tuyến tính tổng quát, máy tăng cường gradient, phân tích Bayes ngây thơ, tích chập mô hình và xgboost. Các thuật toán không giám sát bao gồm mô hình phân cấp tổng quát, mô hình k-means và PCA. Ngoài ra còn có Word2vec cho phân tích văn bản. Phiên bản ổn định mới nhất còn có AutoML – Học máy tự động, một tính năng rất thú vị mà chúng ta sẽ khám phá trong bài viết này.
  3. Chức năng hỗ trợ học máy: Phân tích hiệu suất và tìm kiếm lưới siêu tham số.
  4. Sản phẩm hóa, MapReduce và Đám mây: Hỗ trợ sản phẩm hóa trong môi trường Java; triển khai cụm bằng Hadoop/Spark (Sparkling Water); triển khai trong môi trường đám mây (Azure, AWS, Databricks, v.v.).

Trong bài viết này, chúng ta sẽ sử dụng h2o làm một thuật toán nâng cao cho học máy chuỗi thời gian. Chúng ta sẽ làm việc với h2o trên môi trường cục bộ và phát triển một mô hình chuỗi thời gian có độ chính xác cao trên tập dữ liệu beer_sales_tbl đã được sử dụng trong các hướng dẫn trước về timetk và sweep. Đây là một bài toán hồi quy học máy giám sát.

Chuẩn bị môi trường và dữ liệu

Chúng ta cần ba gói R:

  • h2o: Thư viện cho các thuật toán học máy.
  • tidyquant: Để lấy dữ liệu và tải các công cụ thuộc hệ sinh thái tidyverse.
  • timetk: Thư viện chuyên dụng cho chuỗi thời gian trong R.

Cài đặt H2O

Khuyến nghị cài đặt phiên bản ổn định mới nhất của h2o trên môi trường Ubuntu.

Tải các gói cần thiết

# Tải các thư viện
library(h2o)        # Thư viện ML tuyệt vời
library(timetk)     # Bộ công cụ làm việc với chuỗi thời gian trong R
library(tidyquant)  # Tải tidyverse và các gói tài chính, dùng để lấy dữ liệu

Lấy và Khám phá Dữ liệu

Chúng ta sẽ sử dụng hàm tq_get() từ gói tidyquant để lấy dữ liệu bán hàng bia, rượu vang và rượu mạnh từ FRED.

# Lấy dữ liệu bán hàng bia, rượu vang, rượu mạnh (tính bằng triệu USD)
beer_sales_tbl <- tq_get(
    "S4248SM144NCEN",
    get = "economic.data",
    from = "2010-01-01",
    to = "2017-10-27")

beer_sales_tbl

## # A tibble: 92 x 2
##          date price
##        <date> <int>
##  1 2010-01-01  6558
##  2 2010-02-01  7481
##  3 2010-03-01  9475
##  4 2010-04-01  9424
##  5 2010-05-01  9351
##  6 2010-06-01 10552
##  7 2010-07-01  9077
##  8 2010-08-01  9273
##  9 2010-09-01  9420
## 10 2010-10-01  9413
## # ... with 82 more rows

Việc trực quan hóa dữ liệu là rất quan trọng, đặc biệt là trong phân tích và dự báo chuỗi thời gian. Chúng ta cũng cần chia dữ liệu thành các tập huấn luyện, xác thực và kiểm tra.

# Vẽ biểu đồ doanh số bia với các tập huấn luyện, xác thực và kiểm tra
beer_sales_tbl %>%
    ggplot(aes(date, price)) +
    # Khu vực huấn luyện
    annotate(
        "text",
        x = ymd("2012-01-01"), y = 7000,
        color = palette_light()[[1]],
        label = "Khu vực Huấn luyện") +
    # Khu vực xác thực
    geom_rect(
        xmin = as.numeric(ymd("2016-01-01")),
        xmax = as.numeric(ymd("2016-12-31")),
        ymin = 0, ymax = Inf, alpha = 0.02,
        fill = palette_light()[[3]]) +
    annotate(
        "text",
        x = ymd("2016-07-01"), y = 7000,
        color = palette_light()[[1]],
        label = "Khu vực\nXác thực") +
    # Khu vực kiểm tra
    geom_rect(
        xmin = as.numeric(ymd("2017-01-01")),
        xmax = as.numeric(ymd("2017-08-31")),
        ymin = 0, ymax = Inf, alpha = 0.02,
        fill = palette_light()[[4]]) +
    annotate(
        "text",
        x = ymd("2017-05-01"), y = 7000,
        color = palette_light()[[1]],
        label = "Khu vực\nKiểm tra") +
    # Dữ liệu
    geom_line(col = palette_light()[1]) +
    geom_point(col = palette_light()[1]) +
    geom_ma(ma_fun = SMA, n = 12, size = 1) +
    # Thẩm mỹ
    theme_tq() +
    scale_x_date(
        date_breaks = "1 year",
        date_labels = "%Y") +
    labs(title = "Doanh số bia: 2007 - 2017",
         subtitle = "Hiển thị các tập Huấn luyện, Xác thực và Kiểm tra")

Quy trình Học máy Chuỗi thời gian với H2O và Timetk

Quy trình làm việc cho dự án học máy chuỗi thời gian này tương tự như trong bài viết trước về timetk và hồi quy tuyến tính. Tuy nhiên, lần này chúng ta sẽ thay thế hàm lm() bằng h2o.autoML() để đạt được độ chính xác cao hơn.

Bước 0: Kiểm tra Dữ liệu

Để bắt đầu phân tích, chúng ta sẽ sử dụng glimpse() để xem cấu trúc ban đầu của beer_sales_tbl.

# Điểm khởi đầu
beer_sales_tbl %>%
    glimpse()
## Observations: 92
## Variables: 2
## $ date  <date> 2010-01-01, 2010-02-01, 2010-03-01, 2010-04-01, 20...
## $ price <int> 6558, 7481, 9475, 9424, 9351, 10552, 9077, 9273...

Bước 1: Mở rộng Đặc trưng Chuỗi thời gian

Hàm tk_augment_timeseries_signature() mở rộng thông tin dấu thời gian thành một tập hợp các đặc trưng sẵn sàng cho học máy. Điều này thêm các cột chứa thông tin chuỗi thời gian vào khung dữ liệu gốc. Kiểm tra nhanh bằng glimpse() cho thấy bây giờ chúng ta có 30 đặc trưng.

# Mở rộng (thêm các cột vào data frame)
beer_sales_tbl_aug <- beer_sales_tbl %>%
    tk_augment_timeseries_signature()

beer_sales_tbl_aug %>% glimpse()
## Observations: 92
## Variables: 30
## $ date      <date> 2010-01-01, 2010-02-01, 2010-03-01, 2010-04-01...
## $ price     <int> 6558, 7481, 9475, 9424, 9351, 10552, 9077, 9273...
## $ index.num <int> 1262304000, 1264982400, 1267401600, 1270080000,...
## $ diff      <int> NA, 2678400, 2419200, 2678400, 2592000, 2678400...
## $ year      <int> 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,...
## $ year.iso  <int> 2009, 2010, 2010, 2010, 2010, 2010, 2010, 2010,...
## $ half      <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,...
## $ quarter   <int> 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2,...
## $ month     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3,...
## $ month.xts <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, ...
## $ month.lbl <ord> January, February, March, April, May, June, Jul...
## $ day       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ hour      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ minute    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ second    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ hour12    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ am.pm     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ wday      <int> 6, 2, 2, 5, 7, 3, 5, 1, 4, 6, 2, 4, 7, 3, 3, 6,...
## $ wday.xts  <int> 5, 1, 1, 4, 6, 2, 4, 0, 3, 5, 1, 3, 6, 2, 2, 5,...
## $ wday.lbl  <ord> Friday, Monday, Monday, Thursday, Saturday, Tue...
## $ mday      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ qday      <int> 1, 32, 60, 1, 31, 62, 1, 32, 63, 1, 32, 62, 1, ...
## $ yday      <int> 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 30...
## $ mweek     <int> 5, 6, 5, 5, 5, 6, 5, 5, 5, 5, 6, 5, 5, 6, 5, 5,...
## $ week      <int> 1, 5, 9, 13, 18, 22, 26, 31, 35, 40, 44, 48, 1,...
## $ week.iso  <int> 53, 5, 9, 13, 17, 22, 26, 30, 35, 39, 44, 48, 5...
## $ week2     <int> 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1,...
## $ week3     <int> 1, 2, 0, 1, 0, 1, 2, 1, 2, 1, 2, 0, 1, 2, 0, 1,...
## $ week4     <int> 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 1, 1, 1, 1,...
## $ mday7     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...

Bước 2: Chuẩn bị Dữ liệu cho H2O

Chúng ta cần định dạng lại dữ liệu để h2o có thể đọc được. Đầu tiên, loại bỏ các cột không cần thiết (như cột ngày) và các cột có giá trị thiếu. Sau đó, chuyển đổi các biến thứ tự thành dạng nhân tố thông thường. Các thao tác này có thể thực hiện hiệu quả bằng cách sử dụng gói dplyr.

beer_sales_tbl_clean <- beer_sales_tbl_aug %>%
    select_if(~ !is.Date(.)) %>%
    select_if(~ !any(is.na(.))) %>%
    mutate_if(is.ordered, ~ as.character(.) %>% as.factor)

beer_sales_tbl_clean %>% glimpse()
## Observations: 92
## Variables: 28
## $ price     <int> 6558, 7481, 9475, 9424, 9351, 10552, 9077, 9273...
## $ index.num <int> 1262304000, 1264982400, 1267401600, 1270080000,...
## $ year      <int> 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,...
## $ year.iso  <int> 2009, 2010, 2010, 2010, 2010, 2010, 2010, 2010,...
## $ half      <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,...
## $ quarter   <int> 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 2,...
## $ month     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3,...
## $ month.xts <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, ...
## $ month.lbl <fctr> January, February, March, April, May, June, Ju...
## $ day       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ hour      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ minute    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ second    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ hour12    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ am.pm     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ wday      <int> 6, 2, 2, 5, 7, 3, 5, 1, 4, 6, 2, 4, 7, 3, 3, 6,...
## $ wday.xts  <int> 5, 1, 1, 4, 6, 2, 4, 0, 3, 5, 1, 3, 6, 2, 2, 5,...
## $ wday.lbl  <fctr> Friday, Monday, Monday, Thursday, Saturday, Tu...
## $ mday      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ qday      <int> 1, 32, 60, 1, 31, 62, 1, 32, 63, 1, 32, 62, 1, ...
## $ yday      <int> 1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 30...
## $ mweek     <int> 5, 6, 5, 5, 5, 6, 5, 5, 5, 5, 6, 5, 5, 6, 5, 5,...
## $ week      <int> 1, 5, 9, 13, 18, 22, 26, 31, 35, 40, 44, 48, 1,...
## $ week.iso  <int> 53, 5, 9, 13, 17, 22, 26, 30, 35, 39, 44, 48, 5...
## $ week2     <int> 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1,...
## $ week3     <int> 1, 2, 0, 1, 0, 1, 2, 1, 2, 1, 2, 0, 1, 2, 0, 1,...
## $ week4     <int> 1, 1, 1, 1, 2, 2, 2, 3, 3, 0, 0, 0, 1, 1, 1, 1,...
## $ mday7     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...

Tiếp theo, chúng ta sẽ chia dữ liệu thành các tập huấn luyện, xác thực và kiểm tra dựa trên phạm vi thời gian.

# Chia thành các tập huấn luyện, xác thực và kiểm tra
train_tbl <- beer_sales_tbl_clean %>% filter(year < 2016)
valid_tbl <- beer_sales_tbl_clean %>% filter(year == 2016)
test_tbl  <- beer_sales_tbl_clean %>% filter(year == 2017)

Bước 3: Huấn luyện Mô hình H2O

Đầu tiên, chúng ta khởi tạo h2o để kích hoạt máy ảo Java cần thiết.

h2o.init()        # Khởi động h2o

## Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         46 minutes 4 seconds 
##     H2O cluster version:        3.14.0.3 
##     H2O cluster version age:    1 month and 5 days  
##     H2O cluster name:           H2O_started_from_R_mdancho_pcs046 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   3.51 GB 
##     H2O cluster total cores:    4 
##     H2O cluster allowed cores:  4 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         Algos, AutoML, Core V3, Core V4 
##     R Version:                  R version 3.4.1 (2017-06-30)
h2o.no_progress() # Tắt thanh tiến trình

Chuyển đổi dữ liệu sang định dạng H2OFrame để h2o có thể xử lý.

# Chuyển đổi sang đối tượng H2OFrame
train_h2o <- as.h2o(train_tbl)
valid_h2o <- as.h2o(valid_tbl)
test_h2o  <- as.h2o(test_tbl)

Đặt tên cho biến mục tiêu và các biến dự báo.

# Đặt tên cho h2o
y <- "price"
x <- setdiff(names(train_h2o), y)

Chúng ta sẽ sử dụng h2o.automl() để thử nghiệm các mô hình hồi quy khác nhau trên dữ liệu.

  • x = x: Tên các cột đặc trưng.
  • y = y: Tên cột mục tiêu.
  • training_frame = train_h2o: Tập huấn luyện, bao gồm dữ liệu từ 2010-2016.
  • validation_frame = valid_h2o: Tập xác thực, bao gồm dữ liệu năm 2016, được sử dụng để tránh mô hình bị quá khớp.
  • leaderboard_frame = test_h2o: Mô hình sẽ được xếp hạng dựa trên hiệu suất MAE trên tập kiểm tra.
  • max_runtime_secs = 60: Tham số này giới hạn thời gian tính toán của h2o. Để đảm bảo mô hình chạy được, chúng ta chấp nhận giảm độ chính xác.
  • stopping_metric = "deviance": Sử dụng độ lệch làm chỉ số dừng, điều này có thể cải thiện MAPE của kết quả.
# có thể sử dụng bất kỳ mô hình nào, nhưng ở đây là mô hình hồi quy tuyến tính
automl_models_h2o <- h2o.automl(
    x = x,
    y = y,
    training_frame = train_h2o,
    validation_frame = valid_h2o,
    leaderboard_frame = test_h2o,
    max_runtime_secs = 60,
    stopping_metric = "deviance")

Tiếp theo, chúng ta trích xuất mô hình dẫn đầu (leader model).

# Trích xuất mô hình dẫn đầu
automl_leader <- automl_models_h2o@leader

Bước 4: Dự đoán

Sử dụng hàm h2o.predict() để tạo ra các dự đoán trên tập dữ liệu kiểm tra.

pred_h2o <- h2o.predict(
    automl_leader, newdata = test_h2o)

Bước 5: Đánh giá Hiệu suất

Có nhiều cách để đánh giá hiệu suất mô hình. Ở đây, chúng ta sẽ sử dụng hàm h2o.performance(). Hàm này cung cấp các chỉ số mặc định thường được sử dụng để so sánh các mô hình hồi quy, bao gồm Sai số Căn bậc hai Trung bình (RMSE) và Sai số Tuyệt đối Trung bình (MAE).

h2o.performance(
    automl_leader, newdata = test_h2o)
## H2ORegressionMetrics: gbm
## 
## MSE:  340918.3
## RMSE:  583.8821
## MAE:  467.8388
## RMSLE:  0.04844583
## Mean Residual Deviance :  340918.3

Chỉ số chúng ta quan tâm nhất là Sai số Phần trăm Tuyệt đối Trung bình (MAPE), chỉ số này không có trong kết quả trên. Tuy nhiên, chúng ta có thể dễ dàng tính toán nó. Chúng ta sẽ xem xét sai số trên tập kiểm tra (giá trị thực tế so với giá trị dự đoán).

# Phân tích sai số trên tập kiểm tra
error_tbl <- beer_sales_tbl %>%
    filter(lubridate::year(date) == 2017) %>%
    add_column(
        pred = pred_h2o %>% as.tibble() %>% pull(predict)) %>%
    rename(actual = price) %>%
    mutate(
        error     = actual - pred,
        error_pct = error / actual)

error_tbl
## # A tibble: 8 x 5
##         date actual      pred     error    error_pct
##       <date>  <int>     <dbl>     <dbl>        <dbl>
## 1 2017-01-01   8664  8241.261  422.7386  0.048792541
## 2 2017-02-01  10017  9495.047  521.9534  0.052106763
## 3 2017-03-01  11960 11631.327  328.6726  0.027480989
## 4 2017-04-01  11019 10716.038  302.9619  0.027494498
## 5 2017-05-01  12971 13081.857 -110.8568 -0.008546509
## 6 2017-06-01  14113 12796.170 1316.8296  0.093306142
## 7 2017-07-01  10928 10727.804  200.1962  0.018319563
## 8 2017-08-01  12788 12249.498  538.5016  0.042109915

Để so sánh, chúng ta tính toán một số chỉ số đo lường sai số.

error_tbl %>%
    summarise(
        me   = mean(error),
        rmse = mean(error^2)^0.5,
        mae  = mean(abs(error)),
        mape = mean(abs(error_pct)),
        mpe  = mean(error_pct)) %>%
    glimpse()
## Observations: 1
## Variables: 5
## $ me   <dbl> 440.1246
## $ rmse <dbl> 583.8821
## $ mae  <dbl> 467.8388
## $ mape <dbl> 0.03976961
## $ mpe  <dbl> 0.03763299

Bước 6: Trực quan hóa Kết quả Dự đoán

Cuối cùng, chúng ta trực quan hóa kết quả dự đoán.

beer_sales_tbl %>%
    ggplot(aes(x = date, y = price)) +
    # Dữ liệu - Màu cam
    geom_point(col = palette_light()[1]) +
    geom_line(col = palette_light()[1]) +
    geom_ma(
        n = 12) +
    # Dự đoán - Màu tím
    geom_point(
        aes(y = pred),
        col = palette_light()[2],
        data = error_tbl) +
    geom_line(
        aes(y = pred),
        col = palette_light()[2],
        data = error_tbl) +
    # Thẩm mỹ
    theme_tq() +
    labs(
        title = "Dự báo Doanh số Bia: h2o + timetk",
        subtitle = "H2O đạt độ chính xác cao nhất, MAPE = 3.9%")

Kết quả Cuối cùng

Tỷ lệ MAPE của h2o + timetk tốt hơn so với các phương pháp trong các bài viết trước:

  • timetk + h2o: MAPE = 3.9% (bài viết này)
  • timetk + hồi quy tuyến tính: MAPE = 4.3% (Bộ công cụ phân tích chuỗi thời gian — timetk)
  • sweep + ARIMA: MAPE = 4.3% (Bộ công cụ phân tích chuỗi thời gian — sweep)

Một câu hỏi thú vị có thể đặt ra là: liệu việc lấy trung bình các dự đoán từ ba phương pháp khác nhau có cải thiện độ chính xác hay không?

Lưu ý: Độ chính xác của học máy chuỗi thời gian không phải lúc nào cũng vượt trội so với các kỹ thuật như ARIMA, bao gồm cả các phương pháp được triển khai trong prophet (công cụ dự báo của Facebook) và các phương pháp GARCH. Việc kiểm tra các phương pháp khác nhau và lựa chọn công cụ phù hợp là trách nhiệm của các nhà khoa học dữ liệu.

Thẻ: R H2O Timetk Học Máy chuỗi thời gian

Đăng vào ngày 18 tháng 10 lúc 19:58