Part1 なにはなくともEDA＆ベースモデルつくりから！
Shimpei Ikeno
2022-07-12
本連載の目的：実践的な”多系列”時系列予測モデルの解き方を紹介
みなさんはじめまして。NRIのデータサイエンティスト、時系列予測プラクティスチームの池野です。Wikipediaによれば、時系列とは、“ある現象の時間的な変化を、連続的に（または一定間隔をおいて不連続に）観測して得られた値の系列（一連の値）のこと”をさします。時系列予測の大家であるRob J. Hyndman教授によれば、予測とは、“将来を、過去のデータや影響を与える将来のイベントなどの活用可能な情報に基づいて、できるだけ正確に見通すこと”とあります。したがって、時系列予測は、時間的な変化の観測結果に基づき将来をできるだけ正確に見通す取組といえましょう。時系列予測モデルは、そのような時間的変化の観測結果からパターンを見つけ出しモデル化したものといえます。
ビジネスの世界では、時系列予測というテーマを取り扱う場合、一つの時系列だけを扱えばいいというケースはほとんどありません。会社全体の売上だけでなく、事業別、店舗別、商品別など、様々な単位に分割されたたくさんの系列の予測が必要になります。よくある教科書的な解説では一つ一つの時系列の予測については詳しく書いてあるのですが、たくさんの（階層的な構造がある、あるいはお互いに関係がある可能性がある）時系列をどう取り扱うべきなのかを解説しているソースはほとんどなく、いざビジネス課題に適用しようとすると困ってしまうケースがよくあります。
そこで、本ブログでは、そういったたくさんの時系列（＝多系列）を予測するモデルをどのようにビジネス課題に適用していくのかを連載形式でご紹介していこうと思います。時系列モデルにあまり詳しくないという方も、ある程度知識はあるという方も、１から１０まで丁寧に紹介してまいりますので、ぜひご一読頂ければ幸いです。
時系列予測モデルの現状と課題
さて、時系列予測モデルというとどんな手法があるでしょうか？ARIMAやETSといった古典的な統計モデルから、GBDTなど機械学習を応用したモデル、Facebook（Meta）の開発したProphet、はたまた流行りのDeep Learningまで、多種多様な手法が現在進行形で研究・開発されています。一方で、最新のコンペでも依然としてシンプルな統計モデルのほうが精度が良かったりということもありますし、様々なクラウドベンダや専門企業から提供されているAutoMLの枠組みを使っても思ったより精度が出ないという声もよく聞きます。どんなに素晴らしいモデルでも使い方や使いどころを間違えればその効果は得られません。どんなモデルも簡単に試せるようになった今の時代だからこそ、データサイエンティストは、モデルの特性を正しく理解し、時系列データの特徴を踏まえながら最適なモデルを選定できる目を養う必要があります。
また、説明変数の投入は、一般的なビジネスケースでよく扱われる週次や月次のデータはデータポイントが極端に少なく（例：月次データは単年間12地点、１０年でも１２０地点しかない。）、変数を投入すればするほど精度が上がるということでもないため、丁寧に扱わなければいけない重要な論点です。しかし、説明変数の選択方法について具体的な方法論をアドバイスしているソースはほとんどないように思います。
さらに、単系列の予測ではなく多系列の予測となると問題はさらに複雑化していきます。具体的には、以下のように、たくさんの商品のお店ごとの販売量を予測するようなケースです。
こちらは、“店番＿品番”で構成されているデータになります。もちろん店によっても商品によっても販売傾向はまったく異なりますが、店×品番すべての組み合わせで予測が必要ということになるとその数は簡単に1000を超えますので、例えば1000の時系列に対して5つのモデルを適用しようとすると5000、さらに訓練セットと検証セットに分けるとその2倍、クロスバリデーションしようとすると5倍、とあっという間にモデルの数が増えていきますので、闇雲にいろいろなモデルや変数を試すのは非効率かつ非生産的です。時系列予測のエキスパートになるためには、たくさんのモデルや変数セットを効率的に試行錯誤できるフレームワークを習得する必要があります。
上記の課題認識を踏まえ、本連載では、以下の3点を皆様にお伝えしたいと思います。
- どういうときにどういうモデルを選ぶべきなのか
- 変数作成や選択の考え方
- 効率的なモデリングフレームワーク
本連載の（わかりやすい）目標：（昔の）Kaggleに勝つ！
前述の通り、本ブログではできるだけ実践的な手法をご紹介したいと思っていますが、我々が普段扱うデータはもちろんお客様のデータだったりしますので公開が難しいのが実情です。一方、何か具体的な取組課題がないと読むほうも（書くほうも）なかなかやる気がでませんよね。そこで、今回はKaggleのデータをうまく活用しつつできるだけ我々が実際にプロジェクトで扱うような形でデータ分析のエッセンスをお届けしたいと思います。
対象とするのはこちら(Walmart Recruiting - Store Sales Forecasting | Kaggle)のコンペです。Walmartが2014年に開いた、複数小売店舗の各部門の販売量を予測するというものです。実は、このコンペはリクルーティングを目的としたもので直接的な報酬があるわけでもなく、参加チームも700弱と控え目なのですが、データボリューム的にもハンドルしやすく（数MB程度なのでお手元のラップトップでも十分ハンドルできます。）プロモーションなどの説明変数データもいろいろ入っており、分析・考察の余地も大きいと判断し対象としました。
さて、せっかくコンペに参加するのですから、目標順位は決めたいですよね。2014年といえばまだXGBoostもDeep Learningも時系列ではほとんど活用されてなかったし、ここ数年の技術の進歩は信じられないものがありますよね。ということで、今回の目標は
1位/688チーム
とします！もちろん、過去のコンペなので公開されているノートブックなどでいいスコアを出している例もあるとは思いますが、私たちが確認した限りトップのスコアを越える精度を出しているものはない（はず）なので、こちらを越えられれば少なくとも自力でやったということは皆様にもご理解いただけるかと思います。（とはいえ、あとで思いのほか1位をとることのむずかしさを思い知ることになるのですが、、、それは後々ご紹介できればと思います。）
データ
データは以下の4種類です。
- train.csv: 訓練用のデータです。データの中身はこんな感じ。
## # A tibble: 6 × 5
## Store Dept Date Weekly_Sales IsHoliday
## <dbl> <dbl> <date> <dbl> <lgl>
## 1 1 1 2010-02-05 24924. FALSE
## 2 1 1 2010-02-12 46039. TRUE
## 3 1 1 2010-02-19 41596. FALSE
## 4 1 1 2010-02-26 19404. FALSE
## 5 1 1 2010-03-05 21828. FALSE
## 6 1 1 2010-03-12 21043. FALSE- test.csv: テスト用のデータです。カラムはtrainと同じで、Weekly_Sales だけがないデータです。
- stores.csv: お店の情報が入ったデータです。店単位でマッチする情報になります。
## # A tibble: 6 × 3
## Store Type Size
## <dbl> <chr> <dbl>
## 1 1 A 151315
## 2 2 A 202307
## 3 3 B 37392
## 4 4 A 205863
## 5 5 B 34875
## 6 6 A 202505- features.csv: お店×日付単位でキャンペーン情報などいろいろな情報が入っています。
## # A tibble: 6 × 12
## Store Date Temperature Fuel_Price MarkDown1 MarkDown2 MarkDown3
## <dbl> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2010-02-05 42.3 2.57 NA NA NA
## 2 1 2010-02-12 38.5 2.55 NA NA NA
## 3 1 2010-02-19 39.9 2.51 NA NA NA
## 4 1 2010-02-26 46.6 2.56 NA NA NA
## 5 1 2010-03-05 46.5 2.62 NA NA NA
## 6 1 2010-03-12 57.8 2.67 NA NA NA
## # … with 5 more variables: MarkDown4 <dbl>, MarkDown5 <dbl>, CPI <dbl>,
## # Unemployment <dbl>, IsHoliday <lgl>最後に、サブミッションファイルの形式についても確認しましょう。
## # A tibble: 6 × 2
## Id Weekly_Sales
## <chr> <dbl>
## 1 1_1_2012-11-02 0
## 2 1_1_2012-11-09 0
## 3 1_1_2012-11-16 0
## 4 1_1_2012-11-23 0
## 5 1_1_2012-11-30 0
## 6 1_1_2012-12-07 0このように、“Store” “Dept” “Date”がアンダーバーでつながったものを一つのIdとして、売上を予測するというのが最終的に必要なデータ形式になります。
評価指標
評価指標はWMAEとなっています。基本的にはMAE（Mean Absolute Error）ですが、重みづけ(Weights)があります。
細かい定義はリンクを確認いただくとして、“Holiday Week”の時だけエラーが5倍されるという点が重要です。クリスマスやサンクスギビングなどは非常に多くの売上が期待され、そこでのエラーは他の時期に比べてビジネス上の重要性が高いため、それをモデルに反映したいという意図だと思います。こういったビジネスの意図を反映したカスタムメトリックの設定は我々が実際にモデルを組んでいく際も考慮すべき事柄ですね。
EDA
さて、データと評価指標の紹介が終われば次にくるのはEDA（探索的データ分析）ですね。実際には今回ご紹介すること以上にもっと様々な分析を実施したのですが、すべて紹介していると膨大になってしまうので、今回はまずTrain/Test データに着目し、特にポイントとなる部分だけかいつまんでご紹介したいと思います。
Train データについて
まずはTrainデータの確認です。R PackageであるSkimrを使ってTrainデータをSkimしてみましょう。
| Name | Piped data |
| Number of rows | 421570 |
| Number of columns | 6 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| factor | 3 |
| logical | 1 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Date | 0 | 1 | 2010-02-05 | 2012-10-26 | 2011-06-17 | 143 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Store | 0 | 1 | FALSE | 45 | 13: 10474, 10: 10315, 4: 10272, 1: 10244 |
| Dept | 0 | 1 | FALSE | 81 | 1: 6435, 2: 6435, 3: 6435, 4: 6435 |
| Id | 0 | 1 | FALSE | 3331 | 1_1: 143, 1_1: 143, 1_1: 143, 1_1: 143 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| IsHoliday | 0 | 1 | 0.07 | FAL: 391909, TRU: 29661 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Weekly_Sales | 0 | 1 | 15981.26 | 22711.18 | -4988.94 | 2079.65 | 7612.03 | 20205.85 | 693099.4 | ▇▁▁▁▁ |
Dateは2010年2月～12年の10月までのデータなので、最大で約2年半分の訓練用のデータがあるということですね。StoreとDeptについてはそれぞれ45のお店と81のDeptが対象となります。データとしてはStoreとDeptで分けて週ごとの販売を予測することになりますので、いったんその単位をIdとして確認すると、全部で3331のIdがあるということがわかります。
続いてIsHolidayを見てみると、約7％が該当するということのようです。ちなみに、こちらについてはより具体的に対象とする休暇期間についてKaggleに記載があります。
- Super Bowl: 12-Feb-10, 11-Feb-11, 10-Feb-12, 8-Feb-13
- Labor Day: 10-Sep-10, 9-Sep-11, 7-Sep-12, 6-Sep-13
- Thanksgiving: 26-Nov-10, 25-Nov-11, 23-Nov-12, 29-Nov-13
- Christmas: 31-Dec-10, 30-Dec-11, 28-Dec-12, 27-Dec-13
この4つですね。詳しくは後ほどご紹介しますが、同じ休暇であってもChristmasとSuper BowlではSalesに与える影響がまったく異なりますので、それぞれの休暇は分けて考える必要があります。
最後にWeekly_Salesですが、よく見てみるとどうやらマイナスの売上というのもありそうです。スーパーで売上がマイナスになるというのは考えづらいですが、これ、実際の顧客データを扱っていると本当によくあります。返品とか、もともと予定で入力してしまった分の戻しとか、いろいろな理由でマイナスになっていたりということもありますので注意が必要です。数値の範囲としては大体～2万くらいの範囲に収まっていますが、一部非常に大きい数値（最大約70万）もあるという感じで、右側に長い分布になっていることが想像できますね。
さて、引き続きTrainデータについてもう少しだけ見てみましょう。予測単位であるIdごとにデータがどれだけきちんとそろっているかみてみます。
## # A tibble: 4 × 3
## ndates n prop
## <fct> <int> <dbl>
## 1 (0,40] 3835 0.009
## 2 (40,80] 4800 0.011
## 3 (80,142] 32555 0.077
## 4 (142,143] 380380 0.902143データポイントが完全なデータなのですが、それは全体の約9割にとどまり、残りの10％は期間的な欠損があるということがわかります。こちらもビジネス視点で考えると、期間の途中で新しいお店が開いたり、新しいカテゴリができたり、カテゴリの定義が変更になったりということはよくありますので、データセットとして実践的なものになっていると思い餡巣。先述の通り時系列データはデータポイントが少なくなりがちですので、このような欠損をどう処理するかというのも一つの論点ですね。
Test データについて
Testデータについても同様に概要をSkimしていきましょう。
| Name | Piped data |
| Number of rows | 115064 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| factor | 3 |
| logical | 1 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Date | 0 | 1 | 2012-11-02 | 2013-07-26 | 2013-03-15 | 39 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Store | 0 | 1 | FALSE | 45 | 13: 2836, 4: 2803, 19: 2799, 2: 2797 |
| Dept | 0 | 1 | FALSE | 81 | 1: 1755, 2: 1755, 3: 1755, 4: 1755 |
| Id | 0 | 1 | FALSE | 3169 | 1_1: 39, 1_1: 39, 1_1: 39, 1_1: 39 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| IsHoliday | 0 | 1 | 0.08 | FAL: 106136, TRU: 8928 |
基本的にはTrainデータと構造は変わらないようですが、Dateについては2012年の11月2日の13年7月26日まで、39地点の予測が必要になります。当たり前ですがWeekly_Sales は予測対象なので存在しないですね。Idについては3169とやや少なくなっています。
TrainデータのIdが3331でTestデータのIdが3169ですから、TrainのほうがId数が多いことがわかります。Trainのほうが多いこと自体にはTrainの一部が実際の予測には使われないだけで問題はないはずですが、一応念のためTrainとTestを統合してIdの数を確認しておきましょう。
| Name | Piped data |
| Number of rows | 536634 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| factor | 1 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| Id | 0 | 1 | FALSE | 3342 | 1_1: 182, 1_1: 182, 1_1: 182, 1_1: 182 |
出力するとIdの数は3342（>3331）に増えていることがわかります。ということは、TestにしかないIdが11個は含まれていることになりますね？TrainにないデータをTestで予測するのは時系列モデルではふつう不可能ですから、Id単位に時系列モデルを作っても絶対にサブミットできませんね。。。どうにかしてその辺の問題は解消しなければいけませんが、そちらについて詳しくはモデル構築のところで改めてご紹介しましょう。
Padding & imputation
さて、これまでのところで、TrainデータにはId単位で見ると部分的に欠損があったり、TestにはあってTrainにはないIdがあったりするということがわかりました。したがって、そのまま通常の時系列モデルを適用しようとすると様々なエラーがかえってくることが予見できます。ベースモデル作成に向けて、ここではいったん便宜上簡単に抜けているデータをNAで埋めて（Padding）それらしい値で補完（Imputation）します。
まず、Trainデータで日付的にレコードが欠損しているものをNAで埋めてしまいましょう。以降、時系列の処理には基本的にmodeltime(The Tidymodels Extension for Time Series Modeling • modeltime)というパッケージを利用します。modeltimeは時系列の前処理からtidymodelsの枠組を活用したモデリングまでスムーズに行うことができるパッケージで、とても使いやすいのでおすすめです。データ分析といえばPythonという風潮が強まっていますが、特に時系列分析であればRob J Hyndman教授の一連のパッケージや今回使用するmodeltimeなど非常に使いやすいツールがそろっていますので、ぜひ試してもらえればと思います。でもDeep Learning系のモデルは使えないでしょ？なんて声が聞こえてきそうですが、、、、そちらは後ほどご紹介します。
full_tr <- tr %>%
# create Id
mutate(Id = paste(Store, Dept,sep = "_")) %>%
# padding
group_by(Id) %>%
select(Id, Date, Weekly_Sales) %>%
pad_by_time(
Date,
.start_date = min(tr$Date),
.end_date = max(tr$Date)
) %>%
ungroup() %>%
# change negative to 0
mutate(Weekly_Sales = ifelse(Weekly_Sales<0, 0, Weekly_Sales))ここからは少しCodeの紹介もしていきたいと思います。trはtrainデータで、trに対してまずStoreとDeptをくっつけたIdを作り、Id単位でTrainデータ全体の最初と最後の日付で不足しているDateがある場合NAで埋めるという処理をしています。ついでにWeekly_Salesがマイナスのものを0にしました。
| Name | Piped data |
| Number of rows | 476333 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| Date | 1 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Id | 0 | 1 | 3 | 5 | 0 | 3331 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Date | 0 | 1 | 2010-02-05 | 2012-10-26 | 2011-06-17 | 143 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Weekly_Sales | 54763 | 0.89 | 15981.47 | 22711.03 | 0 | 2079.65 | 7612.03 | 20205.85 | 693099.4 | ▇▁▁▁▁ |
n_missingが０なので、どのIdも欠損のないTrainデータを作ることができました。さきほど見たようにTestだけに存在するIdもあるはずですが、Trainデータに何の情報もないので、Submissionファイルを作る時に0で埋めてしまうのが簡単かと思います。
ベースモデル（Snaive)
欠損のない完全なTrainデータができたところで、いよいよモデリングに入っていきます。ここでは、ベースモデルとしてSnaiveモデルを作成しましょう。SnaiveとはSeasonal naiveモデルのことをさしまして、今回のケースでいえば単純に前年同月同週の数値を予測値として使うというモデルになります。前期の値を予測として使うというのはビジネス上もいまだに本当によく使われておりますので、ベースラインとしてはふさわしいでしょう。
ここでもし対象が一つの時系列であれば単純にモデルを適用していくだけなので苦労はないですが、今回は3000以上の時系列に対して一度にその処理を行いたいというのがミソになります。その場合、大きな方針として
- 一つ一つの時系列に対してモデル適用を繰り返す（Iterative）
- 時系列全体に対して一つのモデルを適用する（Global）
の2つがあります。今回は１．になりますので、それに合わせたデータセット作りから始めます。
# set forecast horizon
horizon <- 39
# create nested data
nested_data_tbl <- tr_imputed %>%
group_by(Id) %>%
extend_timeseries(
.id_var = Id,
.date_var = Date,
.length_future = horizon
) %>%
nest_timeseries(
.id_var = Id,
.length_future = horizon
) %>%
split_nested_timeseries(
.length_test = horizon
)細かい説明はmodeltimeを参照して頂くこととして、上記処理でIdごとの時系列がネストされた状態で保存されたデータを作成することができます。こちらに対してモデルを適用していくのですが、modeltimeはparsnipという、こちらもデータサイエンス界のレジェンド、Max Kuhn博士が開発した素晴らしいフレームワークを活用しています。いろんなモデルを同じSyntaxで適用可能にするというのがおおまかなコンセプトでして、以下のようにmodelとrecipeを別々に定義してworkflowで統合してfitしていくという構造になっています。はじめは扱いが難しいかもしれませんが、慣れると非常に使いやすく、様々なモデル作成・評価を一度に行えたりして非常に便利ですので、modeltimeのフレームワークに則って作成していきましょう。
# set model spec
model_snaive <- naive_reg(
seasonal_period = 52
) %>%
set_engine("snaive")
# create recipe
rec_snaive <- recipe(Weekly_Sales ~ Date,
extract_nested_train_split(nested_data_tbl))
# create workflow
wflw_snaive <- workflow() %>%
add_model(model_snaive) %>%
add_recipe(rec_snaive)まず、modelとrecipeを定義しましょう。modelはparsnipの中からnaive_reg()を選択して、今回は週次データですのでseasonal_period は52でOKです。recipeに特段の工夫はありませんが、iterative modelの場合、使うデータをちょっと特殊な指定の仕方をしなければいけません(extract_nested_train_split)。上記2つをworkflowで統合すれば準備完了です。たくさんの時系列を一遍に処理するとかなり時間がかかることがありますので、まず少量のサンプルでモデルが想定通り回るかどうか確認してから全体に進みましょう。まずは特定のIdだけ抜き出します。
set.seed(235)
id_sample <- sample(nested_data_tbl$Id %>% unique(), size = 6)サンプリングされたデータについてモデルをfitします。
try_sample_tbl <- nested_data_tbl %>%
filter(Id %in% id_sample) %>%
modeltime_nested_fit(
model_list = list(
wflw_snaive
),
control = control_nested_fit(
verbose = TRUE,
allow_par = TRUE
)
)エラーがないかどうか確認し、精度評価やグラフで問題がないか確認しましょう。こういった操作が本の数行のコードで様々なモデルに適用できるというのがmodeltimeの便利なところですね。
try_sample_tbl %>% extract_nested_error_report()## # A tibble: 0 × 4
## # … with 4 variables: Id <chr>, .model_id <int>, .model_desc <chr>,
## # .error_desc <chr>try_sample_tbl %>%
extract_nested_test_accuracy() %>%
table_modeltime_accuracy()| Id | .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
|---|---|---|---|---|---|---|---|---|---|
| 26_25 | 1 | SNAIVE | Test | 893.76 | 13.55 | 0.86 | 14.24 | 1288.86 | 0.35 |
| 29_34 | 1 | SNAIVE | Test | 2185.42 | 22.83 | 1.17 | 20.29 | 2696.7 | 0.32 |
| 44_25 | 1 | SNAIVE | Test | 19.24 | 55.11 | 0.84 | 52.89 | 30.86 | 0.09 |
| 45_67 | 1 | SNAIVE | Test | 1052.62 | 14.77 | 0.57 | 14.1 | 2019.59 | 0.58 |
| 6_1 | 1 | SNAIVE | Test | 4931.23 | 20.43 | 1.1 | 17.28 | 9950.62 | 0.11 |
| 8_19 | 1 | SNAIVE | Test | 631.38 | 122.37 | 5.13 | 60.99 | 1381.9 | 0 |
try_sample_tbl %>%
extract_nested_test_forecast() %>%
group_by(Id) %>%
plot_modeltime_forecast(.interactive = FALSE, .facet_ncol = 3)
エラーは一つもないようですし、グラフを見る限り指示通り52週前のデータを適用しているようです。では、全体に対してモデルをfitしていきます。結構時間がかかるのでparallel processingを適用しましょう。
parallel_start(6)
nested_modeltime_tbl <- nested_data_tbl %>%
modeltime_nested_fit(
model_list = list(
wflw_snaive
),
control = control_nested_fit(
verbose = TRUE,
allow_par = TRUE
)
)6コア/8コアを使用しまして、私のPCで5分ちょっとで終了しました。では、まずエラーとMAEのチェックをしていきます。
# * Review Any Errors ----
nested_modeltime_tbl %>% extract_nested_error_report()## # A tibble: 0 × 4
## # … with 4 variables: Id <chr>, .model_id <int>, .model_desc <chr>,
## # .error_desc <chr>nested_modeltime_tbl %>%
extract_nested_test_accuracy() %>%
head() %>%
table_modeltime_accuracy() | Id | .model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
|---|---|---|---|---|---|---|---|---|---|
| 1_1 | 1 | SNAIVE | Test | 4050.18 | 17.94 | 0.89 | 14.7 | 8977.06 | 0.26 |
| 1_10 | 1 | SNAIVE | Test | 3009.71 | 8.84 | 1.05 | 9.08 | 3988.16 | 0.19 |
| 1_11 | 1 | SNAIVE | Test | 4020.55 | 15.31 | 1.03 | 15.22 | 4959.1 | 0.47 |
| 1_12 | 1 | SNAIVE | Test | 1356.33 | 12.64 | 1.13 | 12.32 | 1675.36 | 0.12 |
| 1_13 | 1 | SNAIVE | Test | 2607.15 | 6.33 | 1.06 | 6.62 | 3113.97 | 0.37 |
| 1_14 | 1 | SNAIVE | Test | 1742.82 | 11.7 | 1.09 | 12.27 | 2232.6 | 0.34 |
エラーはなさそうですね。ＭＡＥは大きく外れているものもありますがまあベースモデルなので気にせずいきましょう。では、次にTrainデータ全部を使ってrefitしていきます。
nested_best_refit_tbl <- nested_modeltime_tbl %>%
modeltime_nested_refit(
control = control_refit(
verbose = TRUE,
allow_par = TRUE
)
)大体3分ちょっとで終了しました。では先ほどと同じようにプロットしてみましょう。
set.seed(234)
id_sample <- sample(nested_data_tbl$Id %>% unique(), size = 6)
nested_best_refit_tbl %>%
extract_nested_future_forecast() %>%
filter(Id %in% id_sample) %>%
group_by(Id) %>%
plot_modeltime_forecast(.facet_ncol = 3)
これだけでもそれなりにいい予測をしているようにも見えますが、例えば3_2や34_35のように、全体としてトレンドが見られるときには少し物足りない予測になっています。
ベースモデルの提出、スコアの確認
では、最後にサブミッション用にファイルを作っていきます。
results <- nested_best_refit_tbl %>% extract_nested_future_forecast()
results <- results %>% mutate(Id = paste(Id, .index,sep = "_")) %>%
select(Id, .value)
results %>% head()## # A tibble: 6 × 2
## Id .value
## <chr> <dbl>
## 1 1_1_2010-02-05 24924.
## 2 1_1_2010-02-12 46039.
## 3 1_1_2010-02-19 41596.
## 4 1_1_2010-02-26 19404.
## 5 1_1_2010-03-05 21828.
## 6 1_1_2010-03-12 21043.submission %>%
left_join(results, by = "Id") %>%
mutate(Weekly_Sales = case_when(!is.na(.value) ~ .value,
is.na(.value) ~ 0
)
) %>% select(Id, Weekly_Sales) %>% write_csv("submission_snaive.csv")今回は予測できていなかったところは0で埋めてしまいました。ファイルを提出しまして、結果を確認すると、、、、
170/688 チーム
Private score: 3026.42451
Public score: 2944.19927
という結果となりました。単純に52期前のデータをあてはめているだけですが、全体の上位約4分の1ということで、意外と悪くない結果ですね。2014年のコンペということで時系列の処理方法が十分浸透してなかったり、Id単位に絞るとデータポイントの制限がかなりあるので機械学習系モデルの適用が難しかったり、原因はいろいろあるかもしれません。Snaiveは基本中の基本ですので、まずはこちらをベースラインとしてどこまで上を目指せるか、引き続き頑張っていこうと思います。
次回は、機械学習系モデルの適用を検討していこうと思いますが、その前段として非常に重要な特徴量について考察していきたいと思います。次回更新にご期待ください！！
