時間序列分析工具箱——sweep

来源:https://www.cnblogs.com/xuruilong100/archive/2018/07/29/9383972.html
-Advertisement-
Play Games

[TOC] 翻譯自《Demo Week: Tidy Forecasting with sweep》 原文鏈接:www.business science.io/code tools/2017/10/25/demo_week_sweep.html 時間序列分析工具箱——sweep 的用途 正如 包之於 ...


目錄

翻譯自《Demo Week: Tidy Forecasting with sweep》

原文鏈接:www.business-science.io/code-tools/2017/10/25/demo_week_sweep.html

時間序列分析工具箱——sweep

sweep 的用途

正如 broom 包之於 stats 包,sweep 包用來簡化使用 forecast 包的工作流。本教程將逐一介紹常用函數 sw_tidysw_glancesw_augmentsw_sweep 的用法。

sweeptimetk 帶來的額外好處是,如果 ts 對象是由 tbl 對象轉換來的,那麼在預測過程中日期和時間信息會以 timetk 索引的形式保留下來。一句話概括:這意味著我們最終可以在預測時使用日期,而不是 ts 類型數據使用的規則間隔數字日期。

載入包

本教程要使用到四個包:

  • sweep:簡化 forecast 包的使用
  • forecast:提供 ARIMA、ETS 和其他流行的預測演算法
  • tidyquant:獲取數據併在後臺載入 tidyverse 系列工具
  • timetk:時間序列數據處理工具,用來將 tbl 轉換成 ts
# Load libraries
library(sweep)      # Broom-style tidiers for the forecast package
library(forecast)   # Forecasting models and predictions package
library(tidyquant)  # Loads tidyverse, financial pkgs, used to get data
library(timetk)     # Functions working with time series

數據

我們使用 timetk 教程中數據——啤酒、紅酒和蒸餾酒銷售數據,用 tidyquant 中的 tq_get() 函數從 FRED 獲取。

# Beer, Wine, Distilled Alcoholic Beverages, in Millions USD
beer_sales_tbl <- tq_get(
    "S4248SM144NCEN",
    get = "economic.data",
    from = "2010-01-01",
    to = "2016-12-31")

beer_sales_tbl
## # A tibble: 84 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 74 more rows

可視化數據是一個好東西,這有助於幫助我們瞭解正在使用的是什麼數據。可視化對於時間序列分析和預測尤為重要。我們將使用 tidyquant 畫圖工具:主要是用 geom_ma(ma_fun = SMA,n = 12) 來添加一個周期為 12 的簡單移動平均線來瞭解趨勢。我們還可以看到似乎同時存在著趨勢性(移動平均線以近似線性的模式增長)和季節性(波峰和波谷傾向於在特定月份發生)。

# Plot Beer Sales
beer_sales_tbl %>%
    ggplot(aes(date, price)) +
    geom_line(col = palette_light()[1]) +
    geom_point(col = palette_light()[1]) +
    geom_ma(ma_fun = SMA, n = 12, size = 1) +
    theme_tq() +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
    labs(title = "Beer Sales: 2007 through 2016")

現在你對我們要分析的時間序列有了直觀的感受,那麼讓我們繼續!

Demo:forecast + sweep 的簡化預測工作流

我們將聯合使用 forecastsweep 來簡化預測分析。

關鍵想法:使用 forecast 包做預測涉及到 ts 對象,用起來並不簡潔。對於 stats 包來說有 broom 來簡化使用;forecast 包就用 sweep

目標:我們將用 ARIMA 模型預測未來 12 個月的數據。

STEP 1:創建 ts 對象

使用 timetk::tk_ts()tbl 轉換成 ts,從之前的教程可以瞭解到這個函數有兩點好處:

  1. 這是一個統一的方法,實現與 ts 對象的相互轉換。
  2. 得到的 ts 對象包含 timetk_idx 屬性,是一個基於初始時間信息的索引。

下麵開始轉換,註意 ts 對象是規則時間序列,所以要設置 startfreq

# Convert from tbl to ts
beer_sales_ts <- tk_ts(
    beer_sales_tbl,
    start = 2010,
    freq = 12)

beer_sales_ts
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct
## 2010  6558  7481  9475  9424  9351 10552  9077  9273  9420  9413
## 2011  6901  8014  9833  9281  9967 11344  9106 10468 10085  9612
## 2012  7486  8641  9709  9423 11342 11274  9845 11163  9532 10754
## 2013  8395  8888 10109 10493 12217 11385 11186 11462 10494 11541
## 2014  8559  9061 10058 10979 11794 11906 10966 10981 10827 11815
## 2015  8398  9061 10720 11105 11505 12903 11866 11223 12023 11986
## 2016  8540 10158 11879 11155 11916 13291 10540 12212 11786 11424
##        Nov   Dec
## 2010  9866 11455
## 2011 10328 11483
## 2012 10953 11922
## 2013 11139 12709
## 2014 10466 13303
## 2015 11510 14190
## 2016 12482 13832

檢查 ts 對象具有 timetk_idx 屬性。

# Check that ts-object has a timetk index
has_timetk_idx(beer_sales_ts)
## [1] TRUE

OK,這對後面要用的 sw_sweep() 很重要。下麵我們就要建立 ARIMA 模型了。

STEP 2A:ARIMA 模型

我們使用 forecast 包里的 auto.arima() 函數為時間序列建模。

# Model using auto.arima
fit_arima <- auto.arima(beer_sales_ts)

fit_arima
## Series: beer_sales_ts 
## ARIMA(3,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##           ar1     ar2     ar3     sar1    drift
##       -0.2498  0.1079  0.6210  -0.2817  32.1157
## s.e.   0.0933  0.0982  0.0925   0.1333   5.8882
## 
## sigma^2 estimated as 175282:  log likelihood=-535.49
## AIC=1082.97   AICc=1084.27   BIC=1096.63

STEP 2B:簡化模型

就像 broom 簡化 stats 包的使用一樣,我麽可以使用 sweep 的函數簡化 ARIMA 模型。下麵介紹三個函數:

  • sw_tidy():用於檢索模型參數
  • sw_glance():用於檢索模型描述和訓練集的精確度度量
  • sw_augment():用於獲得模型殘差

sw_tidy

sw_tidy() 函數以 tibble 對象的形式返回模型參數。

# sw_tidy - Get model coefficients
sw_tidy(fit_arima)
## # A tibble: 5 x 2
##    term   estimate
##   <chr>      <dbl>
## 1   ar1 -0.2497937
## 2   ar2  0.1079269
## 3   ar3  0.6210345
## 4  sar1 -0.2816877
## 5 drift 32.1157478

sw_glance

sw_glance() 函數以 tibble 對象的形式返回訓練集的精確度度量。可以使用 glimpse 函數美化顯示結果。

# sw_glance - Get model description and training set accuracy measures
sw_glance(fit_arima) %>%
    glimpse()
## Observations: 1
## Variables: 12
## $ model.desc <chr> "ARIMA(3,0,0)(1,1,0)[12] with drift"
## $ sigma      <dbl> 418.6665
## $ logLik     <dbl> -535.4873
## $ AIC        <dbl> 1082.975
## $ BIC        <dbl> 1096.635
## $ ME         <dbl> 1.189875
## $ RMSE       <dbl> 373.9091
## $ MAE        <dbl> 271.7068
## $ MPE        <dbl> -0.06716239
## $ MAPE       <dbl> 2.526077
## $ MASE       <dbl> 0.4989005
## $ ACF1       <dbl> 0.02215405

sw_augument

sw_augument() 函數返回的 tibble 表中包含 .actual.fitted.resid 列,有助於在訓練集上評估模型表現。註意,設置 timetk_idx = TRUE 返回初始的日期索引。

# sw_augment - get model residuals
sw_augment(fit_arima, timetk_idx = TRUE)
## # A tibble: 84 x 4
##         index .actual   .fitted    .resid
##        <date>   <dbl>     <dbl>     <dbl>
##  1 2010-01-01    6558  6551.474  6.525878
##  2 2010-02-01    7481  7473.583  7.416765
##  3 2010-03-01    9475  9465.621  9.378648
##  4 2010-04-01    9424  9414.704  9.295526
##  5 2010-05-01    9351  9341.810  9.190414
##  6 2010-06-01   10552 10541.641 10.359293
##  7 2010-07-01    9077  9068.148  8.852178
##  8 2010-08-01    9273  9263.984  9.016063
##  9 2010-09-01    9420  9410.869  9.130943
## 10 2010-10-01    9413  9403.908  9.091831
## # ... with 74 more rows

我們可以可視化訓練數據上的殘差,看一下數據中有沒有遺漏的模式沒有被髮現。

# Plotting residuals
sw_augment(fit_arima, timetk_idx = TRUE) %>%
    ggplot(aes(x = index, y = .resid)) +
    geom_point() + 
    geom_hline(yintercept = 0, color = "red") + 
    labs(title = "Residual diagnostic") +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
    theme_tq()

STEP 3:預測

使用 forecast() 函數做預測。

# Forecast next 12 months
fcast_arima <- forecast(fit_arima, h = 12)

一個問題是,預測結果並不“tidy”。我們需要數據框形式的預測結果,以便應用 tidyverse 的功能,然而預測結果是 forecast 類型的,一種基於 ts 的對象。

class(fcast_arima)
## [1] "forecast"

STEP 4:用 sweep 簡化預測

我們使用 sw_sweep() 簡化預測結果,一個額外的好處是,如果 forecast 對象有 timetk 索引,我們可以用它返回一個日期時間索引,不同於 ts 對象的規則索引。

首先要確認 forecast 對象有 timetk 索引,這需要在使用 sw_sweep() 時設置 timetk_idx 參數。

# Check if object has timetk index 
has_timetk_idx(fcast_arima)
## [1] TRUE

現在,使用 sw_sweep() 來簡化預測結果,它會在內部根據 time_tk 構造一條未來時間序列索引(這一步總是會被執行,因為我們在第 1 步中用 tk_ts() 構造了 ts 對象)註意:這意味著我們最終可以在 forecast 包中使用日期(不同於 ts 對象中的規則索引)!

# sw_sweep - tidies forecast output
fcast_tbl <- sw_sweep(fcast_arima, timetk_idx = TRUE)

fcast_tbl
## # A tibble: 96 x 7
##         index    key price lo.80 lo.95 hi.80 hi.95
##        <date>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 2010-01-01 actual  6558    NA    NA    NA    NA
##  2 2010-02-01 actual  7481    NA    NA    NA    NA
##  3 2010-03-01 actual  9475    NA    NA    NA    NA
##  4 2010-04-01 actual  9424    NA    NA    NA    NA
##  5 2010-05-01 actual  9351    NA    NA    NA    NA
##  6 2010-06-01 actual 10552    NA    NA    NA    NA
##  7 2010-07-01 actual  9077    NA    NA    NA    NA
##  8 2010-08-01 actual  9273    NA    NA    NA    NA
##  9 2010-09-01 actual  9420    NA    NA    NA    NA
## 10 2010-10-01 actual  9413    NA    NA    NA    NA
## # ... with 86 more rows

STEP 5:比較真實值和預測值

我們可以使用 tq_get() 來檢索實際數據。註意,我們沒有用於比較的完整數據,但我們至少可以比較前幾個月的實際值。

actuals_tbl <- tq_get(
    "S4248SM144NCEN",
    get = "economic.data",
    from = "2017-01-01",
    to = "2017-12-31")

註意,預測結果放在 tibble 中,可以方便的實現可視化。

# Visualize the forecast with ggplot
fcast_tbl %>%
    ggplot(
        aes(x = index, y = price, color = key)) +
    # 95% CI
    geom_ribbon(
        aes(ymin = lo.95, ymax = hi.95), 
        fill = "#D5DBFF", color = NA, size = 0) +
    # 80% CI
    geom_ribbon(
        aes(ymin = lo.80, ymax = hi.80, fill = key), 
        fill = "#596DD5", color = NA, 
        size = 0, alpha = 0.8) +
    # Prediction
    geom_line() +
    geom_point() +
    # Actuals
    geom_line(
        aes(x = date, y = price), color = palette_light()[[1]],
        data = actuals_tbl) +
    geom_point(
        aes(x = date, y = price), color = palette_light()[[1]], 
        data = actuals_tbl) +
    # Aesthetics
    labs(
        title = "Beer Sales Forecast: ARIMA", x = "", y = "Thousands of Tons",
        subtitle = "sw_sweep tidies the auto.arima() forecast output") +
    scale_x_date(
        date_breaks = "1 year",
        date_labels = "%Y") +
    scale_color_tq() +
    scale_fill_tq() +
    theme_tq()

我們可以研究測試集上的誤差(真實值 vs 預測值)。

# Investigate test error 
error_tbl <- left_join(
    actuals_tbl,
    fcast_tbl,
    by = c("date" = "index")) %>%
    rename(
        actual = price.x, pred = price.y) %>%
    select(date, actual, pred) %>%
    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  8601.815   62.18469  0.007177365
## 2 2017-02-01  10017 10855.429 -838.42908 -0.083700617
## 3 2017-03-01  11960 11502.214  457.78622  0.038276439
## 4 2017-04-01  11019 11582.600 -563.59962 -0.051147982
## 5 2017-05-01  12971 12566.765  404.23491  0.031164514
## 6 2017-06-01  14113 13263.918  849.08191  0.060163106
## 7 2017-07-01  10928 11507.277 -579.27693 -0.053008504
## 8 2017-08-01  12788 12527.278  260.72219  0.020388035

並且,我們可以做簡單的誤差度量。MAPE 接近 4.3%,比簡單的線性回歸模型略好一點,但是 RMSE 變差了。

# Calculate test error metrics
test_residuals <- error_tbl$error
test_error_pct <- error_tbl$error_pct * 100 # Percentage error

me   <- mean(test_residuals, na.rm=TRUE)
rmse <- mean(test_residuals^2, na.rm=TRUE)^0.5
mae  <- mean(abs(test_residuals), na.rm=TRUE)
mape <- mean(abs(test_error_pct), na.rm=TRUE)
mpe  <- mean(test_error_pct, na.rm=TRUE)

tibble(me, rmse, mae, mape, mpe) %>%
    glimpse()
## Observations: 1
## Variables: 5
## $ me   <dbl> 6.588034
## $ rmse <dbl> 561.4631
## $ mae  <dbl> 501.9144
## $ mape <dbl> 4.312832
## $ mpe  <dbl> -0.3835956

您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • GitHub地址: https://github.com/jeromeetienne/jquery-qrcode ...
  • java的(PO,VO,TO,BO,DAO,POJO)解釋 O/R Mapping 是 Object Relational Mapping(對象關係映射)的縮寫。通俗點講,就是將對象與關係資料庫綁定,用對象來表示關係數據。在O/R Mapping的世界里,有兩個基本的也是重要的東東需要瞭解,即VO, ...
  • 載入靜態文件 在一個網頁中,不僅僅只有一個 html 骨架,還需要 css 樣式文件, js 執行文件以及一些圖片 等。因此在 DTL 中載入靜態文件是一個必須要解決的問題。在 DTL 中,使用 static 標簽來載入 靜態文件。要使用 static 標簽,首先需要 {% load static ...
  • 人應該自信點,因為在某個方面,你無人可取代。做事,做人都要有底線,一件事的底線是什麼,做人的底線是什麼,做事的底線要符合做人的底線。這些事都要清楚。 工作要努力,對你最直接的回饋,就是努力工作所應得的報酬。做人要積極上進,(欲望驅使,興趣驅使,職業規劃,人生態度,生活態度驅使等等) 我今年的計劃,是 ...
  • 在使用 RabbitMQ 的時候,作為消息發送方希望杜絕任何消息丟失或者投遞失敗場景。RabbitMQ 為我們提供了兩個選項用來控制消息的投遞可靠性模式。 rabbitmq 整個消息投遞的路徑為: producer->rabbitmq broker cluster->exchange->que... ...
  • 對於我們開發者來說,設計與實現REST API似乎已經成為了我們的日常生活。API現在已經成為了系統間互通的預設方式。AMAZON就是一個有效的使用API進行系統間溝通的最好的例子。在這篇文章中,我將重點討論如何幫助你設計更好的API並避免一些常見的誤區。 ...
  • 首先點擊File-àNew-àWeb [roject-à在Projcet Name里寫項目名-à點擊finish-à會出來一個框,選擇NO,一個javaweb項目就創建好了。具體請看下圖! 配置伺服器連接: 找到Server-à下麵的Tomcat 7.x點擊右鍵-à點擊Configure Serve ...
  • 想著寫一篇hibernate的博文,於是準備從頭開始,從官網下了最新的穩定版本來做講述。 結果利用hibernate自動建表的時候發生下麵這個問題。 我很納悶,之前用低版本一點的沒有發生這個問題啊。 於是,我把必要文件都拷到之前那個hibernate版本是5.0.7的工程中,結果並沒有發生問題。 所... ...
一周排行
    -Advertisement-
    Play Games
  • 移動開發(一):使用.NET MAUI開發第一個安卓APP 對於工作多年的C#程式員來說,近來想嘗試開發一款安卓APP,考慮了很久最終選擇使用.NET MAUI這個微軟官方的框架來嘗試體驗開發安卓APP,畢竟是使用Visual Studio開發工具,使用起來也比較的順手,結合微軟官方的教程進行了安卓 ...
  • 前言 QuestPDF 是一個開源 .NET 庫,用於生成 PDF 文檔。使用了C# Fluent API方式可簡化開發、減少錯誤並提高工作效率。利用它可以輕鬆生成 PDF 報告、發票、導出文件等。 項目介紹 QuestPDF 是一個革命性的開源 .NET 庫,它徹底改變了我們生成 PDF 文檔的方 ...
  • 項目地址 項目後端地址: https://github.com/ZyPLJ/ZYTteeHole 項目前端頁面地址: ZyPLJ/TreeHoleVue (github.com) https://github.com/ZyPLJ/TreeHoleVue 目前項目測試訪問地址: http://tree ...
  • 話不多說,直接開乾 一.下載 1.官方鏈接下載: https://www.microsoft.com/zh-cn/sql-server/sql-server-downloads 2.在下載目錄中找到下麵這個小的安裝包 SQL2022-SSEI-Dev.exe,運行開始下載SQL server; 二. ...
  • 前言 隨著物聯網(IoT)技術的迅猛發展,MQTT(消息隊列遙測傳輸)協議憑藉其輕量級和高效性,已成為眾多物聯網應用的首選通信標準。 MQTTnet 作為一個高性能的 .NET 開源庫,為 .NET 平臺上的 MQTT 客戶端與伺服器開發提供了強大的支持。 本文將全面介紹 MQTTnet 的核心功能 ...
  • Serilog支持多種接收器用於日誌存儲,增強器用於添加屬性,LogContext管理動態屬性,支持多種輸出格式包括純文本、JSON及ExpressionTemplate。還提供了自定義格式化選項,適用於不同需求。 ...
  • 目錄簡介獲取 HTML 文檔解析 HTML 文檔測試參考文章 簡介 動態內容網站使用 JavaScript 腳本動態檢索和渲染數據,爬取信息時需要模擬瀏覽器行為,否則獲取到的源碼基本是空的。 本文使用的爬取步驟如下: 使用 Selenium 獲取渲染後的 HTML 文檔 使用 HtmlAgility ...
  • 1.前言 什麼是熱更新 游戲或者軟體更新時,無需重新下載客戶端進行安裝,而是在應用程式啟動的情況下,在內部進行資源或者代碼更新 Unity目前常用熱更新解決方案 HybridCLR,Xlua,ILRuntime等 Unity目前常用資源管理解決方案 AssetBundles,Addressable, ...
  • 本文章主要是在C# ASP.NET Core Web API框架實現向手機發送驗證碼簡訊功能。這裡我選擇是一個互億無線簡訊驗證碼平臺,其實像阿裡雲,騰訊雲上面也可以。 首先我們先去 互億無線 https://www.ihuyi.com/api/sms.html 去註冊一個賬號 註冊完成賬號後,它會送 ...
  • 通過以下方式可以高效,並保證數據同步的可靠性 1.API設計 使用RESTful設計,確保API端點明確,並使用適當的HTTP方法(如POST用於創建,PUT用於更新)。 設計清晰的請求和響應模型,以確保客戶端能夠理解預期格式。 2.數據驗證 在伺服器端進行嚴格的數據驗證,確保接收到的數據符合預期格 ...