不確実性を考慮した時系列データ予測

Tweet

はじめに

カブクで機械学習エンジニアをしている大串正矢です。今回は不確実性を考慮した時系列データ予測について書きます。

背景

深層学習では点推定と呼ばれる、ある1点の推定は得意です。例えば、ある画像が与えられたときに、それが90%の確率で猫の画像であるというような推定が可能です。しかし、この確率のブレ幅がどの程度かまでは出してくれません。

85〜95%の90%なのか70〜95%の90％なのかでその信頼度は変わってきます。過去のケースなど実際のデータがある場合は確認は可能ですが未来や未知の予測の場合はブレ幅も分かっていることでどの程度信頼しても良いのか判断できます。このようなブレ幅を不確実性と呼びます。

ベイズの考えを取り入れればこの不確実性を取り入れることが可能です。
実は深層学習にこのベイズの考えを取り入れるための手法としてベイズとDropoutが近似できることを示した論文が出ています。

Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning

今回はこの性質を用いて不確実性をどのように表すか見てみます。

ベイズについて

Dropoutはベイズの近似であることを理解するためにまず、ベイズについて理解しておく必要があります。入力データが\(X\)（例えば画像データなど）、予測すべきラベルが\(Y\)（猫などのラベル）、学習によって得られる重みパラメータが\(W\)（学習で得られる重みパラメータ）とすると学習の際のベイズの式は下記のように記述できます。この式で\(p(W)\)は事前分布になります。

\begin{align}
p(W|X,Y) = \frac{p(Y|X,W)p(W)}{\int{p(Y|X,W)p(W)}dW}
\end{align}

上式の分母の項目は全てのデータを用いているため計算コストが高く、深層学習で有効なミニバッチ処理などが適用できません。そこでミニバッチ処理を適用可能な形に変更していきます。
\(W\)がハイパーパラメータ\(\Theta\)から導出されるような事前分布\(q(W|\Theta)\)を考えます。

\begin{align}
q(W|\Theta) = \arg \min_{\Theta}D(q(W|\Theta), p(W|X,Y))
\end{align}

上の式はパラメータ\(\Theta\)から重みパラメータを生成する分布\(q(W|\Theta)\)と学習データ\(X\)と学習ラベル\(Y\)から重みパラメータ\(W\)が生成される分布\(p(W|X,Y)\)の距離が小さくなるようなパラメータ\(\Theta\)を導出する意味になります。
書き換えると

\begin{align}
D(q(W|\Theta), p(W|X,Y)) = KL(q(W|\Theta)||p(W|X,Y))
\end{align}

分布間のカルバック・ライブラー・ダイバージェンスを小さくすることになります。つまり分布\(q(W|\Theta)\)と分布\(p(W|X,Y)\)がどれだけ似ているかを表しています。この距離が最小化されると分布同士が近い性質になっていることになります。

ではこのカルバック・ライブラー・ダイバージェンスを展開して形を変えます。

\begin{align}
KL(q(W|\Theta)||p(W|X,Y)) = \int_{w}{q(W|\Theta) \log\frac{q(W|\Theta)}{p(W|X,Y)}} \
= – \int_{w} q(W|\Theta) \log\frac{p(W|X,Y)}{q(W|\Theta)} \
\end{align}

ここで\(p(W|X,Y)\)だけに着目し式展開して整理します。
ベイズの定理と条件付き確率の定義を利用します。

まずベイズの定理を利用

\begin{align}
p(W|X,Y) = \frac{p(X,Y|W)p(W)}{p(X, Y)} \
= \frac{p(X,Y|W)p(W)}{p(X, Y)}\frac{p(W)}{p(W)} \
\end{align}

条件付き確率の定義を利用

\begin{align}
= \frac{p(X,Y,W)p(W)}{p(W)p(X, Y)}\frac{p(W)}{p(W)} \
= \frac{p(X,Y,W)p(W)}{p(X, Y)}\frac{1}{p(W)} \
= \frac{p(X,Y,W)}{p(X, Y)} \
= \frac{p(X,Y,W)}{p(X, Y)}\frac{p(X, W)}{p(X, W)} \
\end{align}

条件付き確率の定義を利用

\begin{align}
= \frac{p(Y|X,W)p(X, W)}{p(X, Y)} \
\end{align}

\(X\)と\(W\)が独立、\(X\)と\(Y\)が独立と仮定

\begin{align}
= \frac{p(Y|X,W)p(X)p(W)}{p(X)p(Y)} \
= \frac{p(Y|X,W)p(W)}{p(Y)}
\end{align}

ここで求めた式を用いてカルバック・ライブラー・ダイバージェンスに関する式を書き直します。対数を逆数にしたので\(\Theta\)を最大化するように式を変形します。\(p(Y)\)は\(W\)に依存しないため無視します。

\begin{align}
KL(q(W|\Theta)||p(W|X,Y)) = \arg \max_{\Theta}{- \int_{w} q(W|\Theta) \log\frac{p(Y|X,W)p(W)}{q(W|\Theta)}dW} \
\end{align}

上式を分解します。

\begin{align}
\arg \max_{\Theta}{\int_{w} q(W|\Theta) \log p(Y|X,W)dW – \int_{w} q(W|\Theta) \log\frac{p(W)}{q(W|\Theta)}dW} \
\end{align}
\begin{align}
= \arg \max_{\Theta}{\int_{w} q(W|\Theta) \log p(Y|X,W) dW – KL(q(W|\Theta)||p(W))}
\end{align}

上式をEvidence Lower Bound（略称: ELBO）と呼びます。この式によってパラメータを最適化します。

左の項がデータに対して最適化する項です。
\begin{align}
\int_{w} q(W|\Theta) \log p(Y|X,W)dW
\end{align}

右の項がデータに対して正則化する項です。

\begin{align}
KL(q(W|\Theta)||p(W))
\end{align}

ここで重要なのはベイズの式によって事前分布を与えられること、深層学習においてミニバッチが適用可能な形にするためELBOを用いることになります。

ベイズとDropoutの近似について

Dropoutは下記の図のように考えてもらうとネットワークにノイズを付与していることになります。
ベイズ的に考えるとノイズを加えるような事前分布を考慮していることになります。

重みWは下記のような式になります。

\begin{align}
w_{ij} = \theta_{ij}\epsilon_{ij}, \epsilon_{ij} \sim N(\epsilon_{ij}|1, \alpha)
\end{align}

このノイズを考慮したパラメータの最適化は下記のようになります。stochgradは確率的勾配法を表しています。

\begin{align}
{\rm stochgrad}_{\theta} \log p(Y|X, W)
\end{align}

\begin{align}
= {\rm stochgrad}_{\theta} \log p(Y|X, \Theta\hat{\epsilon})
\end{align}

\begin{align}
\hat{\epsilon} \sim N(\epsilon|1, \alpha I)
\end{align}

ベイズによって重みを導出する事前分布をノイズを含むような正規分布として深層学習のモデルを下記のような式で考えると

\begin{align}
{\rm stochgrad}_{\theta} \int N(W|\Theta, \alpha\Theta^2) \log p(Y|X,W) dW
\end{align}

\begin{align}
N(W|\Theta, \alpha\Theta^2) = \prod N(w_{ij}|\theta_{ij}, \alpha\theta^2_{ij})
\end{align}

ここで\({\rm stochgrad}_{\theta} \log p(Y|X, \Theta\hat{\epsilon})\)と上の式が等価であることが表せれば、ドロップアウトはノイズを加える正規分布と等価として扱えることになります。

これを導出する前にreparameterization-trickについて説明します。
現在は重み\(W\)が分布により導出されている形になります。

これだと重み\(W\)が決定的にならず誤差逆伝搬ができません。そこで仲介するようなパラメータを導入してそのパラメータが分布から生成されると仮定してパラメータ\(W\)は決定的にして誤差逆伝搬を可能にすることをreparameterization-trickと言います。

重み\(W\)が得られる分布をreparameterization-trickによってDropoutによるノイズ付与している正規分布に書き換えています。

\begin{align}
{\rm stochgrad}_{\theta} \int N(W|\Theta, \alpha\Theta^2) \log p(Y|X,W) dW
\end{align}

reparameterization-trickを使用

\begin{align}
= {\rm stochgrad}_{\theta} \int N(\epsilon|1, \alpha) \log p(Y|X, \Theta\epsilon) d\epsilon
\end{align}

\begin{align}
= {\rm stochgrad}_{\theta} \log p(Y|X, \Theta\hat{\epsilon})
\end{align}

\begin{align}
\hat{\epsilon} \sim N(\epsilon|1, \alpha I)
\end{align}

これで最初に定義したDropoutによる効果とベイズで定義した式が等価であることが示せました。
ベイズで定義した式の最適化にはELBOを使用しています。
今回定義した式がELBOと等価かどうかも見てみます。

\begin{align}
\int N(W|\Theta, \alpha\Theta^2) \log p(Y|X,W) dW
\end{align}

ここでノイズの項を置き換えます。

\begin{align}
q(W|\Theta, \alpha) = N(W|\Theta, \alpha\Theta^2) = \prod_{ij} N(w_{ij}|\theta_{ij}, \alpha\theta^2_{ij})
\end{align}

置き換えるとELBOの右のデータ項になります。左の正則化項は\(\alpha\)にのみ依存するようなので今回導出するパラメータとは無関係のため無視できます。参考: Dropout as Bayesian procedure

\begin{align}
\int q(W|\Theta, \alpha) \log p(Y|X,W) dW
\end{align}

これでDrooutとベイズが近似可能であり、パラメータの更新のためのELBOも適用可能なことが分かりました。

Dropoutのベイズ近似を利用した不確実性を考慮した時系列データ予測

Deep and Confident Prediction for Time Series at Uberという論文で時系列データに適用した例があるのでそこで利用されている手法を用いて不確実性を考慮した時系列データ予測を行います。

この論文ではMCdropoutという手法を用いて用いてモデルの不確実性を表しています。
MCDropoutを参考にベイズを考慮したシンプルなアルゴリズムを適用します。

• 1: Dropoutの系列を用意。[0.1, 0.2, 0.3..]など
• 2: Dropoutの系列ごとにモデルを学習
• 3: 学習した各モデルの予測値の平均を導出
• 4: 学習した各モデルの予測値と真の値の誤差の平均を導出
• 5: 3を正規分布の平均に、4で導出した誤差を正規分布の分散に適用

時系列データの取得、前処理、学習のコードはこのブログのコードを使用します。

学習モデルの定義部分です。モデルにdropoutの比率を設定できるようにしています。

def create_model(input_dim,
                 time_steps,
                 latent_dim,
                 # データが一つしかないので1しか選べない
                 batch_size=1,
                 model_option='lstm',
                 optimizer='adam',
                 drop_out=0.5,
                ):
    x = Input(shape=(time_steps, input_dim,))

    if model_option == 'lstm':
        h = LSTM(latent_dim, stateful=False, return_sequences=True, dropout=drop_out)(x)
    elif model_option == 'gru':
        h = GRU(latent_dim, stateful=False, return_sequences=True, dropout=drop_out)(x)

    out = Dense(input_dim)(h)

    model = Model(x, out)
    model.summary()

    model.compile(optimizer=optimizer, loss='mean_squared_error', metrics=['mse'])

    return model

Dropoutのリストを作成して各モデルに適用し、予測のリストと実測と予測の誤差の2乗平均のリストを作成しています。

drop_out_list = [0.01, 0.02, 0.03, 0.04]
predict_list = []
var_list = []

for drop_out in drop_out_list:
    model = create_model(input_dim, 
                         time_steps=time_steps,
                         latent_dim=120,
                         model_option='lstm',
                         drop_out=drop_out,
                        )
    model.fit(x, x, epochs=200)
    window = time_steps
    x_test, scaler = prepare_data(X_smooth_test, time_steps)
    predict_test, x_scale_train = predict_model_show_graph(X_test_day_smooth[window + 1:], x_test, scaler, model)
    predict_list.append(predict_test)
    var_list.append(np.average(np.subtract(x_test, predict_test) ** 2))

全体のコードは下記になります。

https://github.com/SnowMasaya/time_series_anomaly_detect_hands_on/blob/master/advanced/time_series_anomaly_detect_keras_uncertainly.ipynb

評価指標

ここで計測すべき指標は下記になります。

• RMSE: 正常な時系列を正確に再現できているかを表す指標です。低いほど性能が良いことを表します。
• Variance（分散）: 予測波形のブレ幅を表します。小さい方がベターですが一概には言えなく変化が大きい部分では予測が難しくなるのでブレ幅が大きくなります。低いほど一般的な性能が良いことを表します。
• Coverage: 予測波形と分散が実波形をどの程度カバーしているかを表します。高いほど性能が良いことを表します。

下記の図はの際に予測値が取りうる範囲を黄色と緑色で示しています。青が実測値を示しています。

Dropout 1,2,3,4%

Dropoutが10,20,30,40%

今回はDropoutが1,2,3,4%のノイズが少ない事前分布を想定した場合と10,20,30,40%のノイズが中程度の事前分布を比較しました。

下記の表は先ほど示した計測指標のテストデータに対する結果を記述しています。Dropoutを大きくするとRMSEが上昇しますがカバレッジは向上しています。ノイズを多く導入することによりデータに対するカバー率が上昇するのは直感的な感覚に近い現象になっています。

最後に

弊社ではソフトウェアと機械学習を用いて不確実な未来を予測できる方を絶賛採用中なので是非、弊社へ応募してください。

参考

https://www.coursera.org/learn/bayesian-methods-in-machine-learning

https://stats.stackexchange.com/questions/199605/how-does-the-reparameterization-trick-for-vaes-work-and-why-is-it-important