NeurIPS2018読みメモ: Deep State Space Models for Time Series Forecasting(その1)

以下の論文を読みます。

Syama Sundar Rangapuram, Matthias Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, Tim Januschowski. Deep State Space Models for Time Series Forecasting. In Proceedings of NIPS 2018. https://papers.nips.cc/paper/8004-deep-state-space-models-for-time-series-forecasting

※ キャラクターの原作とは関係ありません。解釈誤りは筆者に帰属します。問題点がありましたらご指摘ください。

f:id:cookie-box:20180305231302p:plain:w60

この論文のタイトルは「時系列予測のためのディープ状態空間モデル」だって。ちなみに NeurIPS2018 にはこの論文とは別に「条件なし単語生成のためのディープ状態空間モデル( Deep State Space Models for Unconditional Word Generation )」っていうタイトルめっちゃかぶってる論文もあるから注意な。

f:id:cookie-box:20180305232608p:plain:w60

時系列予測? そんなのどんな時系列かによるでしょう? 検証に利用しているデータは何ですか?

f:id:cookie-box:20180305231302p:plain:w60

なんで俺が説教されるの!? えっと、この論文では定性評価のための検証と定量評価のための検証とやってるみたいだけど、定性評価パートでは人工データのパラメータを復元できるかを検証してるみたいだから、現実の時系列データをつかってるのは定量評価パートだな。定量評価パートでつかってるのは、

  • "electricity": これは370人の顧客の時間ごとの電力消費量?
  • "traffic": サンフランシスコの高速道路の963の車線の時間ごとの占有率? 混み具合ってこと?
どっちも1日の中での「季節性」があるってかいてある。夜の方が電気がよく消費されるとか、通勤時間は道路が混んでるとか、そんな感じかな。public なデータみたいだから、ちょっとみてみるか。どこで手に入るんだ? まず以下の論文への参照があるな。なんかこの論文もこの論文で内容気になるけど、この論文が "electricity" と "traffic" データを検証につかってるみたい。で何々… "The data sets electricity and traffic are obtained from the UCI repository" ってあるから UCI リポジトリってところから取ってくるのか。じゃあ最初の論文でそうかいてくれよ…。とりあえず "UCI repository" で検索すれば出てきそうだな。
  • https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014
    • "electricity" の方はたぶんこれだな…時刻のカラムを除くとカラム数が 370 ある。370人の顧客の電力消費量って370人の顧客の電力消費量の合計じゃなくて個々の顧客の電力消費量って意味だったのか…。数値は15分ごとの電力消費量(kW)で、時刻はポルトガル時間で毎年3月と10月にそれぞれ23時間しかない日と25時間ある日があるって。サマータイムだな。
  • https://archive.ics.uci.edu/ml/datasets/PEMS-SF
    • "traffic" の方はこれか。これは 10 分ごとの 963 個のセンサーの測定値みたいだ。PEMS_train ファイルは行ごとに 963 × 144 の行列になってるから、これが1日分のデータ(963個のセンサーの10分ごとの測定値=963×6×24=963×144)で、PEMS_trainlabels ファイルは PEMS_train ファイルの各行が「何曜日のデータだったか」を示してるな(1~7の各数字が月~日曜日を意味する)。曜日を識別する分類問題を意識したデータみたいだ。データ内には深夜にセンサーが止まって数値が欠損している日も2日含まれてるらしい…。
これらのデータの一部をプロットしてみるとこうかな。顧客370人分は多すぎるから5人分だけで、センサー963個は無理だから3個だけで、期間もほんの一部に絞ってあるけど。

f:id:cookie-box:20190121200857p:plain
f:id:cookie-box:20190121210707p:plain
f:id:cookie-box:20180305231302p:plain:w60

"electricity" は人が寝静まってる深夜の時間帯が落ち込んでる感じなのかな? みんな寝てるからパソコンしないとか? "traffic" はよくわかんないけど、やっぱり深夜のが道路が空いてる? 日曜日と比べて水曜日は朝方に道路が混んでるな。でも夜は若干日曜のが混んでる? って考えるとなんか特徴はありそうだけど、全部の曜日のパターンを見抜くのは難しそうだな…。ただこの論文では曜日の分類はしてなくてあくまで数値の予測をしようとしてるんだと思うけど。

f:id:cookie-box:20180305232608p:plain:w60

どちらも10~15分毎に観測されるデータで、異なる顧客/車線から同種のデータが大量に得られるんですね。

f:id:cookie-box:20180305231302p:plain:w60

アブストラクトによると、時系列予測のために状態空間モデルと深層学習を組み合わせたみたいなんだけど、時系列ごとの線形状態空間モデルを RNN で結合してパラメタライズすることで??、状態空間モデルと深層学習のいいとこどりをしてるって…状態空間モデルから受け継いだいいところは、解釈性と data efficiency?

f:id:cookie-box:20180305232608p:plain:w60

そのまま読めばデータを効率よくつかえるということなんでしょうが…カルマンフィルタのゲインのようなイメージでしょうか。カルマンゲインによって新しく観測された情報を必要なだけ取り入れるような。

f:id:cookie-box:20180305231302p:plain:w60

まあ詳しくは後で出てくるだろ。深層学習の方のいいところは複雑なパターンも学習できると。それで提案手法は、学習データが少ししかないときにデータの構造を想定するのにも、学習データが大量にあるときに複雑なパターンを学習するのにも効果的らしい…。Conclusions をみてもかいてあることはあまり変わらないな。提案手法と比較したのは最近の RNN ベースの手法と、matrix factorization 手法と、伝統的な手法だとはかいてある。あと提案手法の利点は季節構造を明示的にモデリングできるとこらしい。…あのさ、この論文ちょいちょい Appendix を参照ってあるんだけど、Appendix なくない?

f:id:cookie-box:20180305232608p:plain:w60

確かにないですね…あ、Proceedings のページの "PDF" ではなくて "Supplemental" というリンク先のファイルには Appendix が付いていますよ。

f:id:cookie-box:20180305231302p:plain:w60

マジで? よかった、メールで問い合わせないといけないかと思った。だったら俺英作文できないからハルナに頼まなきゃだし。

f:id:cookie-box:20180305232608p:plain:w60

春名さんは英語を話せるけど書けないんじゃ…というか春名さんにまず日本語で趣旨が伝えられるでしょうか。「論文の Appendix が見当たらないのですが」という質問をしたいんですよね?

f:id:cookie-box:20180305231302p:plain:w60

駄目だその日本語ハルナに通じないな…。まあ Appendix あったからいいや。Intoroduction を読むと時系列予測でも特に需要予測は現実で重要だってかいてるな。確かに "electricity" と "traffic" も需要予測っぽいな。道路の需要っていうとちょっと違和感あるけど。とにかくそういう予測では正確で最新のデータに基づいた予測値がほしいと。それで、時系列のトレンドや季節性を捉えるモデルの最たる例が状態空間モデルだけど、状態空間モデルは時系列の構造があらかじめ結構わかっているときには有利で、解釈もしやすいし、データ効率?のいい学習もできるけど、一方でじゅうぶんな長さの学習データが必要だし、あと巨大で複雑な時系列データを扱うには説明変数を決めるのが大変すぎだって。さらに、伝統的な状態空間モデルでは時系列をばらばらにフィッティングするから類似の時系列が共有するパターン?を推論できない? そうなの?

f:id:cookie-box:20180305232608p:plain:w60

状態空間モデル自体はもちろん多変量時系列を扱うことができますし、その多変量の分布をトラッキングできますが…「類似の時系列」というのがポイントなんでしょうか。「類似の時系列」というのは明らかに「いろいろな顧客の消費電力の時間変化のパターン」「いろいろな道路の混み具合の時間変化のパターン」のようなものを指していますよね。このような時系列を予測するとき、「人が活動しない夜中は消費電力が少ないだろう」などという特徴を考慮してシステムモデルをつくるんだと思います。でもそのシステムモデルを370人の顧客それぞれに適用したのでは、370人分のデータがある意味があまりないんだと思います。もちろん「顧客Aと顧客Bにの消費電力には強い正の相関がある」などというのは結果的にトラッキングされますけど、「顧客Aと顧客Bにの消費電力が強い正の相関をもつのはどのようなときである」というのを積極的に学習するわけではない、というようなイメージでしょうか…。

f:id:cookie-box:20180305231302p:plain:w60

ふーん? それで逆に深層ニューラルネットワークは高次元の特徴を抽出することができて、あまり人が手間をかけることなく時系列を横断する複雑なパターンを特定できるって。本当かなー。ただこの手のモデルはあまり時系列の構造に何か想定をおいてるものじゃないから、正確な学習データを得るにはやっぱり長期間分の学習データが必要で、解釈もしにくいし、平滑化みたいなことをしたくても織り込みにくいって。

f:id:cookie-box:20180305232608p:plain:w60

「伝統的なモデルは解釈性と可塑性に優れるが手間を要し表現力に欠け、ニューラルネットは手間がかからず表現力に長けるが解釈性と可塑性に欠ける」という導入に幾度とない既視感がある気がしますね。

f:id:cookie-box:20180305231302p:plain:w60

それで便利なモデルを考えたなら別にいいだろ! でこの論文はもちろん状態空間モデルとニューラルネット、とりわけ RNN みたいだけど、を組み合わせるんだけど、「線形状態空間モデルを RNN でパラメタライズする」らしい…これだけじゃ意味わかんないな。でも、個別の時系列は状態空間モデルでモデリングして、時系列間を横断する構造を RNN でモデリングするっぽい?

f:id:cookie-box:20180305232608p:plain:w60

人工データで検証しているということですから、どのような人工データを生成しているかで著者らがどのような構造の学習を意図しているのかわかるでしょうね。

f:id:cookie-box:20180305231302p:plain:w60

提案モデルは状態空間モデルをつかってるから解釈もしやすいし、状態空間モデルに入れるモデルの想定によって少ないデータにも大量のデータにも適応できる。ってデータが少なくても多くてもいいって状態空間モデルパート任せだったのか…。まあだから過学習を防げるって。データが少ないと過学習しやすいから、データが少ないなりの想定を状態空間モデルに入れとけってことだよな…。とまあイントロダクションとしてはこんな感じで、2節は先行研究だけど。

f:id:cookie-box:20180305232608p:plain:w60

短いので確認しておきましょうか。Hyndman et al. や Durbin and Koopman が状態空間モデルについてまとめており、最近では機械学習分野で状態空間モデルに関する研究もみられると。状態空間モデルと RNN の組合せというのはこれまでにも提案されていて、以下の片方もしくは両方に該当するということです。

  • ガウシアンではない複雑な尤度関数に拡張する。
  • MLP によって時間発展に非線形性をもたせる、または状態空間モデルと RNN によって特定した時間変化する遷移行列を組み合わせる。
というか状態空間モデルと深層学習の組合せの例がかなり挙げられていますが、論文上の記述からのみではアイデアがわかりづらいのでタイトルだけ列挙しておきます。
  • Deep Markov Model (DMM) [Krishnan2017, Krishnan2015]
  • Stochastic RNNs [Fraccaro2016]
  • Variational RNNs [Chung2015]
  • Latent LSTM Allocation [Zaheer2017]
  • State-Space LSTM [Zheng2017]
  • 状態空間モデルの教師なし学習 [Karl2017]
  • Kalman Variational Auto-Encoder (KVAE) [Fraccaro2017] ― 本論文での提案手法に最も近い。
  • 上と似たアイデア [Johnson2016]
この下から2番目のモデルが本論文の提案手法と似ているらしいですが、こちらの手法では RNN で固定された K 個のパラメータを学習しているのに対して、本論文の手法では RNN で直接的に状態空間モデルのパラメータを学習するのでハイパーパラメータチューニングが不要とありますが…よくわからないので先に進みましょう。

f:id:cookie-box:20180305231302p:plain:w60

3節は一般的な状態空間モデルによる予測の話らしいけど、表記や前提をつかむために読んでおいた方がいいな。いま  \{ z_{1:T_i}^{(i)}\}_{i=1}^N っていう N 本の時系列があるっぽい。 \{ z_{1:T_1}^{(1)} \}, \; \{ z_{1:T_2}^{(2)} \}, \cdots, \; \{ z_{1:T_N}^{(N)}\} ってことだな。それで  z_t^{(i)} \in \mathbb{R} だな。脚注に、各時系列は等間隔だけど1本目の時系列の時刻 t=1 と2本目の時系列の時刻 t=1 はそろってなくてもいいってあるけど、まあ気にしなくていいかな? それで、それに伴う説明変数もあって  \{ x_{1:T_i + \tau}^{(i)}\}_{i=1}^N こう。こっちの  x_t^{(i)} はベクトルで  x_t^{(i)} \in \mathbb{R}^D だ。なんで  x_t^{(i)} の方が  z_t^{(i)} より \tau だけ長く手元にあるんだって思ったけど、\tau ステップ先が予測対象で、 x_t^{(i)} はその \tau ステップ先までの値が手に入っていて、その状態で  z_{T_i + 1:T_i + \tau}^{(i)} を予測するのか。以下の確率分布を求めたいらしい。

p \left( z_{T_i + 1:T_i + \tau}^{(i)}  \, \middle| \, z_{1:T_i}^{(i)} , \, x_{1:T_i + \tau}^{(i)}   ; \, \Phi \right)
縦棒の右側に他の時系列は入ってこないんだな。でもモデルのパラメータ \Phi は他の時系列とシェアされてるらしい。それで仮定として、  x_{1:T_i}^{(i)}\Phi が given の下では N 本の時系列  \{ z_{1:T_i}^{(i)}\}_{i=1}^N は時系列間で独立とするって。この仮定を置いちゃうと時系列間の相関構造はモデリングできないけど…って、あ!!

f:id:cookie-box:20180305232608p:plain:w60

どうしたんですか?

f:id:cookie-box:20180305231302p:plain:w60

時系列間が独立と仮定しちゃうと相関構造をモデリングはできないけど、パラメータ \Phi が共通だから、時系列間で「統計的な強度」をシェアできないわけじゃない、ってかいてあるんだよ。statistical strength って、PRML の6章の最後で出てきたじゃん!

このときも、ニューラルネットの各出力変数は隠れ層を共有するから互いに「統計的な強度」を借りることができるっていってた。

f:id:cookie-box:20180305232608p:plain:w60

多変量を予測するモデルは、個別にモデリングするよりもまとめてモデリングして共通因子をもたせる方が個々の変数をより確かに予測できるという意味なんでしょうね。

f:id:cookie-box:20180305231302p:plain:w60

それで、状態空間モデルは典型的には1本1本の時系列に対して適用するって。こういう電力消費量の予測みたいな文脈ではってこと? だって気象予報の状態空間モデルでは全然そんなことないよな…。まあともかく、こういうシステムモデルを考える。

 l_t = F_t l_{t-1} + g_t \varepsilon_t , \qquad \varepsilon_t \sim N(0,1)
l_t が状態で F_t が遷移行列で、ノイズもガウシアンだし、カルマンフィルタが適用できる状況だな。それで観測モデルがこう。やっぱり線形でガウシアンな観測を考えてる。
 z_t = y_t + \sigma_t \epsilon_t, \quad y_t = a_t^{\top} l_{t-1} + b_t , \quad \epsilon_t \sim N(0,1)
それで、特定しなきゃいけないパラメータは  \Theta_t = (\mu_0, \Sigma_0, F_t, g_t, a_t, b_t, \sigma_t) で、これを特定するには例えば周辺尤度を最大化するようなパラメータを選ぶ。周辺尤度は以下だな。
 \displaystyle p_{SS}(z_{1:T}| \Theta_{1:T}) \equiv p(z_1|\Theta_1) \prod_{t=2}^{T} p(z_t | z_{1:t-1}, \Theta_{1:t})
それでパラメータが特定できたら z_{1:T} の分布は l_{1:T}積分すれば出るけど、今回の条件では z_{1:T} の分布は l_{1:T}積分しなくてもカルマンフィルタで解析的に求まるな。でも、せっかく N 本の時系列があるのに1本1本予測してたんじゃ時系列間で情報が共有されないから、学習データが少なかったりノイズが大きかったりしても適用できるようにこれから RNN と組み合わせるって。

f:id:cookie-box:20180305232608p:plain:w60

それで4節がその Deep State Space Models…各時刻に  x_t^{(i)} を入力し、各時刻の出力が  \Theta_t^{(i)} となるように RNN を学習するということですか。論文タイトルが Deep State Space Model"s" と複数形なのも N 本の各時系列の状態空間モデルを同時に学習するからだったんですね。ネットワークの学習の損失関数は周辺尤度ですね。RNN は時系列からその時刻までの特徴を抽出するのに適したモデルなので、学習対象を状態空間モデルそのものにするというのは理解できる気がしますが…これ、"electricity" と "traffic" における「説明変数(状態)」とは何なんでしょう?

f:id:cookie-box:20180305231302p:plain:w60

4節端折りすぎだろ! この論文のメインパートなのに! 5節の検証パートは Appendix を合わせてみないとわからないなこれ…。これたぶん"electricity" と "traffic" ではどっちも状態空間モデルでも季節性に基づくモデルを採用してて、トラッキングしてる状態は 31 次元で、そのうち 24 次元は1日の中での時間の季節成分、7 次元は曜日の季節成分をあらわしてる。詳細はこっちの論文みろって。この DeepAR は比較対象にしてる手法の1つだな。

次回があればつづく

electricity のプロット
import numpy as np
import pandas as pd
import datetime

df = pd.read_csv('LD2011_2014.txt/LD2011_2014.txt', sep=';', dtype='object', nrows=600)
df.columns.values[0] = 'time' # 時刻カラムだけ列名がないので付ける

# 時刻カラムの処理
df['time'] = df['time'].map(lambda x: datetime.datetime.strptime(x, '%Y-%m-%d %H:%M:%S'))

# 時刻でないカラムの処理 (このデータは小数点がカンマになっている)
for i in range(1, 371):
    df.iloc[:,i] = df.iloc[:,i].map(lambda x: x.replace(',', '.')).astype(np.float64)
    #if df.iloc[:,i].sum() > 0.0:
    #    print(i)

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 14, 6
rcParams['font.size'] = 16

plt.plot(df['time'], df['MT_191'])
plt.plot(df['time'], df['MT_192'])
plt.plot(df['time'], df['MT_193'])
plt.plot(df['time'], df['MT_194'])
plt.plot(df['time'], df['MT_195'])
plt.legend()
plt.xlabel('time')
plt.ylabel('electricity')
plt.show()
traffic のプロット
import numpy as np

x = []
f = open('PEMS-SF/PEMS_train')
line = f.readline()
n = 0
while line:
    rows = line.strip('[]\n').split(';') # ファイルの各行は大カッコでくくられ、行の境目がセミコロンで、行内が空白で区切られた行列
    x_ = []
    for row in rows:
        x_.append(np.array(row.split(' ')).astype(np.float64))
    x.append(x_)
    n = n + 1
    if n == 11:
        break
    line = f.readline()
f.close()

x = np.array(x)
#print(x)
#print(x.shape)

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 14, 6
rcParams['font.size'] = 16

# 0番目が水曜日で10番目が日曜日( PEMS_trainlabels をみる)
plt.plot(np.arange(0, 24, step=1/6), x[0][0], label='Wednesday Sensor0', color='dodgerblue')
plt.plot(np.arange(0, 24, step=1/6), x[0][1], label='Wednesday Sensor1', color='dodgerblue', linestyle='dashed')
plt.plot(np.arange(0, 24, step=1/6), x[0][2], label='Wednesday Sensor2', color='dodgerblue', linestyle='dotted')
plt.plot(np.arange(0, 24, step=1/6), x[10][0], label='Sunday Sensor0', color='darkorange')
plt.plot(np.arange(0, 24, step=1/6), x[10][1], label='Sunday Sensor1', color='darkorange', linestyle='dashed')
plt.plot(np.arange(0, 24, step=1/6), x[10][2], label='Sunday Sensor2', color='darkorange', linestyle='dotted')
plt.legend()
plt.xlabel('hour')
plt.ylabel('occupancy rate')
plt.show()