M5 Forecasting 参加報告 m5 反省会 Jul 19, 2020 Team Kon https://github.com/Y-oHr-N/m5-forecasting

Yu Ohori (a.k.a. Kon) NS Solutions Corporation (Apr 2017 - ) ● Researcher ● Data Science & Infrastructure Technologies ● System Research & Development Center ● Technology Bureau @Y_oHr_N @Y-oHr-N #SemiSupervisedLearning #AnomalyDetection #DataOps

3

Outline Introduction Overview Difficulties Solutions Discussion Conclusion 4

Introduction kaggle とは？ ● 2010 年に設立されたデータ分析者のためのコミュニティ ○ 300 近くのデータ分析コンペを開催 ○ 1 つのデータ分析コンペに数百ないし数千ものユーザが参加 ● 2017 年に Google によって買収 ● 2017 年にユーザ数が 100 万人を突破 5 活躍の場を求む データ分析者 データ活用に悩む 企業または政府 データ，賞金 モデル

Introduction 内資系企業もデータ分析コンペを開催 ● リクルート、世界最大のデータサイエンティストコミュニティ"Kaggle"と日本企業初の共催となるデータ予測コン ペティション「RECRUIT Challenge - Coupon Purchase Prediction」開催を決定 | リクルートホールディングス ● Kaggleは凄かった! 更に簡単な出品を目指して商品の値段推定精度を改善中 | メルカリエンジニアリング 採用時に kaggle の実績を優遇する企業が増加 ● ヤフー、技術者採用に初年度650万円提示　新卒も可: 日本経済新聞 いくつかの企業で業務中の kaggle への参加を推奨 ● 業務時間のKaggle参加を推奨する新制度導入。DeNA内でのAI技術者育成と活用を強化 | 株式会社ディー・エヌ・ エー【DeNA】 ● RistのKaggleチーム創設に伴い、Kaggle Grandmaster 小野寺和樹氏がアドバイザーとして就任決定 | 株式会社 Rist | Rist Inc. 6

Introduction 7 車載カメラに写る車両の位置・姿勢を予測 （自動車） 遺伝子発現パターンから化合物の作用機序を予測 （医薬） Wikipedia のページ毎の日別 PV を予測 （通信） 店舗・商品毎の日別売上を予測 （小売） 画像や説明文からその商品の需要を予測 （広告） 組織サンプルをもとに前立腺癌の進行度を予測 （医療） 顧客が金融商品を契約するか否かを予測 （銀行） 自動車保険加入者が翌年に保険金請求するか否かを予測 （保険） 鋼板の表面にある欠陥を検出 （鋼材）

Outline Introduction Overview ● data ● task ● metric ● timeline Difficulties Solutions Discussion Conclusion 8

Overview M5 Competition とは？ ● MOFC が開催する 5 回目の時系列予測コンペ ● Walmart がデータを提供 ○ 米国のスーパーマーケットチェーン ● Kaggle で 2 つの部門を同時に開催 ○ Accuracy ■ 5,558 teams ■ 7,022 competitors ○ Uncertainty ■ 909 teams ■ 1,103 competitors ● 日本勢が Accuracy の 300 位以下に多数位置 9 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/163410

Overview - data 3 つの表からなるデータセット（size: 429 MB） ● calendar.csv ○ 暦（年，月，日，曜日） ○ 行事（祝日，スポーツ，宗教等） ○ フードスタンプ（SNAP）の支給日 ● sell_prices.csv ○ 約 5 年分の店舗・商品毎の週別価格 ● sales_train_{validation,evaluation}.csv ○ 各種 ID ○ 約 5 年分の店舗・商品毎の日別売上数 2 つの部門で同一のデータを使用 災害情報等の外部データを利用可能 10 店舗名，商品名，部門名は不明

Overview - task M5 Forecasting - Accuracy ● 次の 28 日間の日別売上数を予測 M5 Forecasting - Uncertainty ● 次の 28 日間の日別売上数を 9 つの分位点で予測 ○ u = 0.005, 0.025, 0.165, 0.250, 0.500, 0.750, 0.835, 0.975, 0.995 11 Level id Aggregation level #series 1 all 1 2 state_id 3 3 store_id 10 4 cat_id 3 5 dept_id 7 6 state_id, cat_id 9 7 state_id, dept_id 21 8 store_id, cat_id 30 9 store_id, dept_id 70 10 item_id 3,049 11 item_id, state_id 9,147 12 item_id, store_id 30,490

Overview - metric 12 Root Mean Squared Scaled Error (RMSSE) ● Y_{t}: 時刻 t の実測値 ● \hat{Y}_{t}: 時刻 t の予測値 Weighted RMSSE (WRMSSE) ● 階層時系列を考慮した Accuracy の評価指標 ● w_{i}: 時系列 i の最終 28 日間の合計売上高 Scaled Pinball Loss (SPL) ● Y_{t}: 時刻 t の実測値 ● Q_{t}(u): 時刻 t の u 分位点における予測値 Weighted SPL (WSPL) ● 階層時系列を考慮した Uncertainty の評価指標 ● w_{i}: 時系列 i の最終 28 日間の合計売上高

Overview - metric 13 何故縮尺を調整して評価する必要があるのか？ ● 時系列の縮尺の違いに左右されずに予測精度を評価したい為 何故重み付けして評価する必要があるのか？ ● 売上高の大きい時系列に対する予測精度を優先的に評価したい為

Overview - timeline 14 検証フェーズ ● train: 11/01/29 (d_1) - 16/04/24 (d_1913) ● test: 16/04/25 (d_1914) - 16/05/22 (d_1941) ● weight: 16/03/28 (d_1886) - 16/04/24 (d_1913) ● 結果を public score として算出 評価フェーズ ● train: 11/01/29 (d_1) - 16/05/22 (d_1941) ● test: 16/05/23 (d_1942) - 16/06/19 (d_1969) ● weight: 16/04/25 (d_1914) - 16/05/22 (d_1941) ● 結果を private score として算出 ● 最終 submission を 1 つ選択 20/03/01 競技開始 20/06/30 競技終了 20/06/01 追加データ公開

Outline Introduction Overview Difficulties ● zero-inflated data ● nonstationary time series ● hierarchical time series ● multi-step time series forecasting Solutions Discussion Conclusion 15

Difficulties 今回の問題で難しい要素は？ ● 零過剰データ ○ 過剰な 0 に引っ張られて予測が低めに偏る ● 非定常時系列 ○ 季節変動，トレンド，イベントの有無のせいで対象期間によって予測精度が大きく変化する ● 階層時系列 ○ 42840 つの時系列を精確に予測する必要がある ● マルチステップ予測 ○ 後段のステップは利用できる情報が限られるため，ラグ特徴の作成で不利を強いられる 以上の要素に対処するための幾つかのアプローチを Discussion から紹介 16

Difficulties - zero-inflated data 店舗・商品単位で見ると売上数の 67 % が 0 に分布 ● ただし，未発売期間の売上数は 0 0 が不自然に 長く続く期間が存在（右図右下） ● 供給障害による在庫切れの可能性 ○ 食糧品をメキシコから調達 ○ 非食糧品を中国から調達 以下の方法で零過剰データに対処 ● tweedie 回帰 ● Croston 法 ● 分類器の併用 17 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/138085

Difficulties - zero-inflated data Croston 法 [Croston, '72] ● 最も古い断続的需要予測手法 ● 新たに需要が発生しない限り，予測値は一定 ● Y_{t}: 時刻 t の実測値 ● \hat{Y}_{t}: 時刻 t の予測値 ● q: 最後に需要が発生した時点からの経過時間 ● \alpha: 平滑化パラメータ 18 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/133644

Difficulties - zero-inflated data TSB 法 [Teunter+, '11] ● Croston 法の派生手法 ● 新たに需要が発生しない限り，予測値は減少 ● Y_{t}: 時刻 t の実測値 ● \hat{Y}_{t}: 時刻 t の予測値 ● \alpha, \beta: 平滑化パラメータ 19 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/133644

Difficulties - nonstationary time series 決定木に基づいたモデルは外挿が不得手（右図） ● 非定常時系列の予測に不向き 以下の方法で非定常時系列に対処 ● 統計モデルを用いてトレンドを除去 ● 時系列を差分系列に変換 20 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/155950 train test

Difficulties - hierarchical time series 21 方法名 内容 モデル作成に使用する時系列 集約による情報の損失 bottom-up 法 最下位の時系列を用いて学習・予測 これらを合計して上位層の時系列を作成 不安定 有 top-down 法 最下位の時系列を用いて学習・予測 これを分解して下位層の時系列を作成 安定 無 middle-out 法 任意の層の時系列を用いて学習・予測 bottom-up 法で上位層の時系列を作成 top-down 法で下位層の時系列を作成 https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/163595

Difficulties - multi-step time series forecasting 22 戦略名 内容 ステップ間の依存関係の把握 再帰的予測による誤差の累積 学習時間 direct multi-step strategy モデルをステップ数だけ作成 各ステップを別々に予測 できない ない 長い recursive multi-step strategy モデルを 1 つ作成 各ステップを再帰的に予測 できる 大きい 短い direct-recursive hybrid strategy モデルをステップ数だけ作成 各ステップを再帰的に予測 できる 小さい 長い multiple output strategy 多出力モデルを 1 つ作成 全ステップを一度に予測 できる ない 長い https://www.kaggle.com/c/m5-forecasting-accuracy/discussion/141972

Outline Introduction Overview Difficulties Solutions ● top solutions ● Kon’s solution Discussion Conclusion 23

Solutions - top solutions 24 順位 チーム名 ラグ特徴 モデル 目的関数 magic multiplier 1 YeonJun IN_STU 非再帰的，再帰的 LGBM tweedie 無 2 Matthias 無 LGBM, N-BEATS asymmetric rmse 無 3 mf 再帰的 DeepAR tweedie 無 4 monsaraida 非再帰的 LGBM tweedie 無 5 Alan Lahoud 非再帰的，再帰的 LGBM poisson 有 magic multiplier: 予測値に掛けるだけで score を大きく改善させる乗数

Solutions - top solutions 25 https://twitter.com/spyrosmakrid/status/1283756501581037571 Twitter で レベル別順位を公開

Solutions - Kon’s solution 26 名前 所属 Kaggle の経験 得意領域 主な担当 chimuichimu (Ichimura) NSSOL, FIN 無 EDA FE Kon (Ohori) NSSOL, SYSRDC 有（1 silver） tabular FE modeling (CB, LGBM) ko10 (Shimizu) NSSOL, STEEL 有 image, text modeling (NN) yamazombie (Yamada) NSSOL, IND 無 metric

Solutions - Kon’s solution 27 利用したサービス・ツール 備考 コミュニケーション Microsoft Teams, Slack 週 1 回の定例以外は Slack を利用 タスク管理 Trello 利用頻度が低く運用に失敗 ナレッジ管理 Scrapbox 日報や議事録，情報共有に利用 実験管理 Google Sheets 競技終盤から利用 コード管理 Git, GitHub データ管理 DVC 競技中盤から利用 環境管理 pip LightGBM の実験再現に苦戦 計算資源 Google Cloud Platform 競技終盤から n1-highmem-16 を利用

28 7 月 5 月 3 月 4 月 6 月 3/9-3/15 チーム結成 ● メンバーを募集 ● 各種サービスを整備 3/15-3/25 問題理解と EDA ● データは？ ● タスクは？ ● 評価指標は？ 3/26-5/2 初 submit ● ベースラインを実装 ● 未経験者を補助 ● Notebooks を調査 ● Discussion を調査 5/3-6/1 特徴・モデル作成 ● アイデアを量産 ● FE ● CB, NN を実装 6/1-6/13 事後分析と CV ● Public LB を事後分析 ● 宗教行事を調査 ● WRMSSE を実装 ● CV を実装 6/14-6/30 実験 ● GCP に移行 ● 人力特徴選択 ● 人力 HPO Solutions - Kon’s solution

Solutions - Kon’s solution 29 M5 Forecasting - Accuracy M5 Forecasting - Uncertainty

Solutions - Kon’s solution 30 前処理 時系列毎に未発売期間を除去 ● 0 の割合を減らす為 ● rolling の計算から除外する為 時系列毎に価格の外れ値を置換 価格を掛けて目的変数を変換 ● WRMSSE に寄せる為 後処理 負の予測値を 0 で置換 価格で割って目的変数を逆変換 特徴生成 9 categorical features ● 5 id features ● 4 calendar features 98 numerical features ● 12 calendar features ● 29 price features ● 57 sales features モデル LGBM ● learning_rate=0.03 ● num_boost_round=875 ● objective="tweedie" recursive multi-step strategy

Solutions - Kon’s solution 31 calendar に関する主な特徴（右図赤） ● event_name_{1,2}, event_type_{1,2} ○ LGBM の内部でエンコーディング ● 8 つの calendar_features ○ attrs = [year, dayofyear, weekofyear, month, quarter, day, weekofmonth, weekday]

Solutions - Kon’s solution 32 price に関する主な特徴（右図緑） ● 20 つの aggregate_features ○ level 2, …, 11 に対して作成 ○ agg_funcs=[mean, std] ● 4 つの expanding_features ○ level 12 に対して作成 ○ agg_funcs=[min, max, mean, std] ● 2 つの pct_change_features ○ periods=[7, 28]

Solutions - Kon’s solution 33 sales に関する主な特徴（右図青） ● 6 つの shift_features ○ level 10, 11, 12 に対して作成 ○ periods=[7, 28] ● 48 つの rollinig_features ○ level 10, 11, 12 に対して作成 ○ windows=[7, 14, 28, 56] ○ agg_funcs=[mean, std] ○ min_periods=1

Solutions - Kon’s solution ハイパパラメータを手動で探索（右図） ● Optuna を参考に探索順・探索空間を決定 ● 以下のように分割して CV を実施 ○ fold 1: valid=[d_1858, …, d_1885] ○ fold 2: valid=[d_1886, …, d_1913] ○ fold 3: valid=[d_1914, …, d_1941] ● fold 毎に early stopping を実施 ○ WRMSSE を監視 ○ best_iteration の平均を再学習時に使用 34 ハイパパラメータ名 探索空間 ベストハイパパラメータ num_leaves 32, 64, …, 2048 512 bagging_fraction 0.70, 0.75, …, 0.90 0.8 bagging_freq 1 1 feature_fraction 0.70, 0.75, …, 0.90 0.8 lambda_l1 0.0 0.0 lambda_l2 0.0001, 0.001, …, 1.0 0.001 min_data_in_leaf 1000, 1500, …, 2500 1500 tweedie_variance_power 1.05, 1.10, …, 1.25 1.2

Outline Introduction Overview Difficulties Solutions Discussion ● ML vs statistics ● magic multiplier Conclustion 35

Discussion - ML vs statistics Makridakis, S., et. al, "Statistical and Machine Learning forecasting methods: Concerns and ways forward," 2018. ● 時系列予測において 8 つの統計モデルと 10 つの機械学習モデルの性能を比較 ● M3 Competition データセットを用いて実験 ● 予測精度と計算時間の両面で統計モデルの方が優位 36 sMAPE

Discussion - ML vs statistics 19 つの統計モデルと 5 つの機械学習モデルを Accuracy のベンチマークとして公開 ● 4 つの統計モデルが銅圏に位置 ○ Exponential Smoothing - bottom-up（ES_bu) ○ Exponential Smoothing with eXplanatory variables (ESX) ○ Average of the two ES methods, the first computed using the top-down approach and the second using the bottom-up approach (Com_tb) ○ AutoRegressive Integrated Moving Average with eXplanatory variables (ARIMAX) 50 位以下の参加者は 1 つか 2 つを除いて全て機械学習モデルを利用 37

Discussion - magic multiplier 何故予測値を 1.04 倍すると public score が大きく改善 [Konstantin Yakovlev, 3/29] したのか？ ● 時系列に上昇トレンドが見られた為 [John Curcio, 4/6] ○ 2016 年第 2 四半期の売上高が前期に比べて増加 [Maxwell, 3/29] ○ 外挿が不得手なモデルを補助 ○ 上昇トレンドは item の増加によるもの，モデルに関係なく score が改善した為，この仮説は不適 ● 検証期間（d_1914 - d_1941）に外れ値が含まれていた為 [miguel perez, 7/3] ○ calendar.csv に記載のない Pentecost (5/15) に予測値が大きく下振れ [Kon, 6/9] ○ 上位層の時系列の帳尻が合うように調整 何故予測値を 0.955 倍すると private score が大きく改善 [Giba, 7/1] したのか？ ● 検証期間に外れ値が含まれていた為 ○ 検証期間のデータをもとに作成したラグ特徴によって予測値が上振れ 38

Outline Introduction Overview Difficulties Solutions Discussion Conclusion 39

Conclusion 業務に向けた教訓は？ ● 受領データのみで予測が満足に行えるのかを検討する ● モデルの良し悪しを評価する際にテストデータを十分に用意する ● モデルの良し悪しにそれほど差がない場合，管理が容易なモデルを選択する ● 分析者・依頼者の双方が解釈できる評価指標を使用する 40

Appendix - articles about LightGBM for beginners 以下の順に閲覧することを推奨 ● 導入 ○ 「初手LightGBM」をする7つの理由 - u++の備忘録 ● xgboost, catboost との比較 ○ PyData Tokyo Meetup #21 LightGBM ● ハイパパラメータ最適化 ○ PyData.Tokyo Meetup #21 講演資料「Optuna ハイパーパラメータ最適化フレームワーク」太田 健 ○ なぜn_estimatorsやepochsをパラメータサーチしてはいけないのか - 天色グラフィティ ● 特徴重要度 ○ 特徴量重要度にバイアスが生じる状況ご存知ですか？ - 学習する天然ニューラルネット ● 学習アルゴリズム ○ Kaggle Masterが勾配ブースティングを解説する - Qiita ● カテゴリカルエンコーディングと欠損値補完 ○ lightgbm カテゴリカル変数と欠損値の扱いについて+α - てばさきさんの自由研究 ○ How to encode categorical features for GBDT - Speaker Deck ○ 機械学習における欠損値補完について考える - rmizutaの日記 ○ 推論時の入力に未知の欠損値がある場合のlightgbmの挙動の確認 - rmizutaの日記 ● その他の文献 ○ GBDTの理解に役立つサイトまとめ - copypasteの日記 41

Appendix - LightGBM tips LightGBM reproducibility | Kaggle ● 使用したコンパイラが異なるとシードを固定しても学習結果を再現できなくなる On impact of num_threads to LGBM training time | Kaggle ● 大規模なマシンを用いて学習する際，全スレッドを使用すると却って学習時間が悪化する LightGBM reproducibility issue | Kaggle ● 列の順番を並び替えると学習結果を再現できなくなる LightGBMError: Do not support special JSON characters in feature name. | Kaggle ● 特徴名にカンマのような特殊文字が含まれていると例外を送出する 42

Appendix - LightGBM tips Three shades of Dark | Kaggle ● bagging を行う際，学習結果が n_threads の値に依存する ● 2.3.2 で修正 Speed-up LightGBM training | Kaggle ● コンパイル時にコンパイラフラグを設定すると学習を高速化できる How to avoid LightGBM memory spike when it starts training | Kaggle ● lgb.Dataset に pd.DataFrame を渡すと construct する際にメモリを大量に消費する does anyone using Lightgbm gpu version? | Kaggle ● GPU 上で学習する際，categorical feature の水準数が極端に大きいと例外を送出する 43