Re:ゼロから始めるML生活

機械学習ド素人が、流行に乗ってニューラルネットとかその他いろいろ勉強してみるブログ

【R】時系列分析の覚書(基本、カルマンフィルタ)

前回はウィナーフィルタまでやって力尽きました。

tsunotsuno.hatenablog.com

そんなわけで、今回もこちらの本をやっていきます。

今回はカルマンフィルタです。 一応位置づけは8章で、基本になっています。

f:id:nogawanogawa:20180715133812j:plain:w500

(これで基本とか、世の中難しすぎだろ…)

カルマンフィルタ

カルマンフィルタとは

何をやっているのかを見失うと訳わかんなくなるので、目的を明確にします。 カルマンフィルタの目的は、

時系列データが与えられたときに、観測値から状態を推定すること

とします。

カルマンフィルタでは、前回同様、線形・ガウス型状態空間モデルを想定します。 式としてはこんな感じです。


x_t = G_t x_{t-1} + w_t


y_t = F_t x_t + v_t

このモデルに従う時系列データで観測値から状態を推定することを考える際にカルマンフィルタを適用することができます。 カルマンフィルタリングのアルゴリズムを書くとこんな感じです。

  1. 時点t-1でのフィルタリング分布:m_{t-1}, C_{t-1}
  2. 時点tでの更新手続き
    • 一期先予測分布
      • (平均) a_t \gets G_tm_{t-1}
      • (分散) R_t \gets G_tC_{t-1}G_{t}^{\mathrm{T}} + W_{t}
    • 一期先予測尤度
      • (平均) f_t \gets F_ta_{t}
      • (分散) Q_t \gets F_tR_{t}F_{t}^{\mathrm{T}} + V_{t}
    • カルマン利得 K_t \gets R_tF_{t}^{\mathrm{T}}Q_t^{-1}
    • 状態の更新
      • (平均) m_t \gets a_t + K_{t} [ y_t - f_t ]
      • (分散) C_t \gets [I-K_tF_t ]  R_t
  3. 時点tでのフィルタリング分布:m_{t}, C_{t}

・・・さっぱりわかりません! ってことで、ちょっと他の方のブログを参考にする反則技を使用します。

qiita.com

カルマンフィルター - ORWiki

非常に参考になりました。 ありがとうございます。

もとになる制御工学なんかで出てくるカルマンフィルタの図がこちら。

f:id:nogawanogawa:20180916164429j:plain:w600

そんでもって、今回のケースに合うようにいろいろ情報を修正した版がこちら。

f:id:nogawanogawa:20180916164509j:plain:w600

これだとなんかわかる感じがします。ちなみに図の上半分は、モデルとしてはあるけれど実際の点線の中身は観測できません。

観測できないxによって決定される観測値yの値を、カルマンフィルタ側(下半分)で推定してあげます。 推定していくと、途中でxの推定値から次のyを推定できるようになります。 このyの推定値と実際のyの観測値を比較して推定の誤差を小さくしていきます。

このときに出てくるのがカルマンゲインと言うやつです。 カルマンゲインでは、本当のxに対してxの推定値との誤差を最も小さくすることを目標に算出します。 すなわちそれは、xの平均二乗誤差を最も小さくなったときに達成され、平均二乗誤差を最小にすることでカルマンゲインが得られます。

ここまで簡単にしてもやっぱりむずいっすね。。。尤度とか分布とか確率モデルを考える前に、普通に実数で考えたほうが良さそうです。

カルマンフィルタを使ってみる

何はともあれなんとなくやってみるのが一番はやいので、やってみます。 この本では、サンプルとしてナイル川の流量のデータを使用した状態の推定をやってみます。

  • 観測値 : y(Nile)
  • 状態 : m

理屈は難しかったんですが、実装は式のとおりに書いていくだけなんで、そこは意外と簡単です。

フィルタリング

フィルタリングなんで、現在までの観測値を使って現在の状態を推定します。

予測

予測なんで、現在までの観測値を使って未来の状態を推定します。 前提として、フィルタリングが完了している(標準とする時刻まで状態の推定が完了している)必要があります。

予測のアルゴリズムを書いておくとこんな感じです。

  1. 時点t+(k-1)でのk-1期先予測分布:a_t(k-1), R_t(k-1)
  2. 時点t+kでの更新手続き
    • k期先予測分布
      • (平均) a_t(k) \gets G_{t+k}a_{k}(k-1)
      • (分散) R_t(k) \gets G_{t+k}R_{t}(k-1)G_{t+k}^{\mathrm{T}} + W_{t+k}
  3. 時点tでのフィルタリング分布:a_{t}(k), R_{t}(k)

平滑化

平滑化なんで、現在までの観測値を使用して過去の状態を推定します。 こちらもフィルタリングがしてある(標準とする時刻まで状態の推定が完了している)前提になります。

平滑化のアルゴリズムとしてはこんな感じです。

  1. 時点t-1でのフィルタリング分布:s_{t+1}, S_{t+1}
  2. 時点tでの更新手続き
    • 平滑化利得 A_t \gets C_tG_{t+1}^{\mathrm{T}}R_{t+1}^{-1}
    • 状態の更新
      • (平均) s_t \gets m_t + A_{t} [ s_{t+1} - a_{t+1} ]
      • (分散) S_t \gets C_t + A_{t}[S_{t+1}-R_{t+1} ]  A_t^{\mathrm{T}}
  3. 時点tでのフィルタリング分布:s_{t}, S_{t}

ローカルレベルモデルでやってみる

なんとな~く状態の推定ができそうだと言うことがわかったところで、今度はちゃんと分析してみます。

対象データ

観測値はナイルを使用していました。図だけ出すとこんな感じです。 データに関する事前の準備については、以前もやったのですがもう一度やってみます。

なんとなく見ていくと、5数要約としてはこんな感じですね。

  • 最小値:456.0
  • 中央値:893.5
  • 平均値:919.4
  • 最大値:1370.0

感覚的には外れ値もなさそうです。 自己相関(acf)としては特定の相関は見られない気がします。 スペクトルを見てみても、相関が見られないような感じですね。

ローカルレベルモデルの定義


x_t = x_{t-1} + w_t


y_t = x_t + v_t

いつもとほとんど同じ式ですね。 GとFが常に1なだけです。 Fが1なので、状態と観測値の間にホワイトノイズしか乗ってきていないことになります。 つまり、ノイズの分の誤差はあれど、それ以外は状態と観測値は一致することを意味しています。

分析

さっきとほとんど同じなんですが、一応分析の全体像を示します。

f:id:nogawanogawa:20180917223221j:plain:w500

上では、著者お手製のカルマンフィルタを使って推定しましたが、今回はRのdlmライブラリを駆使して推定するため、 コードの流れが変わりますが、全体像としてやることは同じです。 ライブラリを使うためのお作法による違いがあるだけです。

Nileのデータがローカルレベルモデルに従うと仮定して、モデルを構築します。 その後、フィルタリング・予測・平滑化によって状態を推定します。

モデルの構築

何はともあれモデルを作ってみます。 この辺はもう書いてあるとおりですね。。。

フィルタリング

ローカルレベルモデルになると若干モデルの支配方程式が変わるので、結果も自ずと変わってきます。 まずはフィルタリングをやってみます。

赤い実線がフィルタリングによる状態の推定の平均、緑の実践が95%信頼区間の境界を示しています。 代替がNileのグラフを覆っているので、そこそこの制度で見積もれているように見受けられます。

予測

次に予測です。

これは、未来の状態を求めているので、まあこうなりますよね。 上下の破線は95%信頼区間を表しています。 時間の経過とともに正確に推定できず、範囲が広がっていく様子がわかります。

平滑化

最後に平滑化です。

一番最後に一つおまけでグラフを付けてみました。 95%信頼区間の推移を色別に示しています。

フィルタリングに比べて平滑化のほうが気持ち範囲が狭まっていることがわかりますでしょうか? 平滑化によって、より詳細に状態が推定されていることがわかります。

結果の吟味

さて、ちゃんと出てきた結果のについて評価します。

観測値yに対して、仮定したモデルが適合しているかどうかを判断することになります。 実際には、様々なモデルや前提に対して状態を推定してそれぞれを相対的に比較することで、現象がどのモデルに従っているかを判断します。

今回はローカルレベルモデルしかやっていないので、基準となる指標を出してみます。

まず尤度について出しています。

次にイノベーション(予測誤差・残差)を出しています。 ここで言うイノベーションeとは観測値yと一期先予測尤度fに対して


e_t = y_t + f_t

で表されます。

これを一期先予測尤度の分散で割ったe_t / \sqrt{Q_t}を規格化イノベーションと言います。

適切にモデルが設定されていれば、規格化イノベーションはノイズの影響による正規分布に従うはずです。 今回はそちらも確認します。

最後に、QQプロットを出しています。 規格化イノベーションeのプロットがy=xの式に近づけば近づくほど正規分布に従っていることを表します。 まあ大体正規分布かなというレベルですね。

そんなわけで、参考値程度ではありますがなんとなく適切にモデルを適用できているように見受けられます。

感想

いやー、難しいです。 次はいろいろ分析をやってみます。

【R】時系列分析の覚書(基本、ウィナーフィルタ)

前回に引き続き、こちらの本をやっていきます。

前回は導入ということで、簡単な分析をなぞってみました。

tsunotsuno.hatenablog.com

今回の範囲はこちらです。

f:id:nogawanogawa:20180715133812j:plain:w500

前回はとりあえず簡単な時系列分析をやってみるといった位置づけでしたが、今回は状態空間モデルの基本をやってみます。

状態空間モデル

時系列データがどんな関係性にあるかを表現する際に、潜在変数があるかないかで大きく二分されます。 潜在変数を使用するモデルの一つが状態空間モデルで、逆に潜在変数を使用せず観測されるデータの関係性を直接表現するモデルの一つとしてARMAモデルが挙げられます。

状態空間モデルを図で示すと下のようになります。

f:id:nogawanogawa:20180721091602j:plain:w500

一方でARMAモデルを図にするとこんなイメージです。

f:id:nogawanogawa:20180908170953j:plain:w400

観測されるデータが何に依存しているかが最大の違いとなっています。

この本では状態空間モデルは次のように説明されています。

状態空間モデルは互いに関連のある系列的なデータを確立的に捉えるモデルの1つです。

確率空間モデルでは、直接観測されるデータ(観測値)と直接には観測されない確率変数(状態)を用います。 観測値や状態には、分析を行う上で解釈が容易になるように量を選ぶ必要があります。 そこで、ざっくり言うと下記のような関係性を仮定します。

状態については、状態は前時点のみに依存する(マルコフ性
観測値については、ある時点の観測値についてはその時点の状態にのみ依存する

これとは対象的に、観測される時系列データだけを見て関係性を求めるモデルをARMAモデルと言ったりします。

状態空間モデルを確率分布で表現すると以下のようになります。

p(x_t | x_{0:t-1}, y_{1:t-1}) = p(x_t | x_{t-1})
p(y_t | x_{0:t}, y_{1:t-1}) = p(y_t | x_{t})

方程式で表すとこんな感じです。

x_t = f(x_{t-1}, w_t) (状態方程式、システム方程式)
y_t = h(x_t, v_t) (観測方程式)

wvのことを、それぞれ状態雑音(システム雑音)観測雑音と言います。

ちなみに、図でも確率分布でも方程式でも意味はすべて同じです

状態の推定

状態空間モデルでの状態を推定することを考えます。

観測値は取得可能ですが、状態は観測不可能です。 そのため、分析に使用する際には、状態を推定することがポイントになります。

フィルタリング(現在の状態の推定)

フィルタリングは過去から現在まで観測された時系列データを使用して、現在の状態を推定します。

  • input : 観測された時系列データ
  • output : 現在の状態

イメージにするとこんな感じです。

f:id:nogawanogawa:20180721143132j:plain:w500

見方を変えるとこんな感じで表現されます。

f:id:nogawanogawa:20180721143109j:plain:w500

これは右方向が潜在変数の時間進行、下方向が観測値の時間進行を表しています。 実際には同じ時間軸に乗っているので、右下の方向に向かって進んでいきます。

フィルタリング分布について細かい数式はおいておいて、結論だけ書くと、

フィルタリング分布: 
p(x_t | y_{1:t}) = p(x_t | y_{1:t-1}) \frac{p(y_t | x_t}{p(y_t | y_{1:t-1})})

一期先予測分布: 
p(x_t | y_{1:t-1}) =
  \int p(x_t | y_{1:t-1}) p(x_{t-1} | y_{1:t-1}) dx_{t-1}

一期先予測尤度: 
p(y_t | y_{1:t-1}) =
  \int p(y_t | x_t) p(x_t | y_{1:t-1}) dx_t

となります。(証明は端折ります。本買って読んでください。)

簡単にまとめると、

フィルタリングでは時刻0の状態から始まり、

  1. 潜在変数が時間進行
  2. 潜在変数の変化に伴い観測値が変化

を繰り返していきます。イメージとしてはこんな感じです。

f:id:nogawanogawa:20180721143149j:plain:w500

このようにフィルタリングでは階段状に、現在の観測値から潜在変数を推定していきます。

予測(未来の状態を推定する)

次に未来の観測値について考えます。 未来なので観測値は手にはいりません。

f:id:nogawanogawa:20180721143209j:plain:w500

この場合は、手に入っている現在の観測値から潜在変数を求めていきます。(証明は端折ります。本買って読んでください。)

予測分布: 
p(x_{t+k} | y_{1:t}) =
  \int p(x_{t+k} | x_{t+k-1}) p(x_{t+k-1} | y_{1:t}) dx_{t+k-1}

イメージはこんな感じで潜在変数だけ時間進行するイメージで潜在変数を推定します。

f:id:nogawanogawa:20180721143221j:plain:w500

平滑化(過去の状態を推定する)

今度は現在より過去の状態を推定することを考えます。

f:id:nogawanogawa:20180721143236j:plain:w500

一旦時刻T(>t)までフィルタリングが完了している状態から考えます。 平滑化分布は下記のように表現されます。(証明は端折ります。本買って読んでください。)

平滑化分布: 
p(x_t | y_{1:T}) = p(x_t | y_{1:t})
\int\frac{p(x_{t+1} | x_t}{p(x_t+1 | y_{1:t})} p(x_{t+1} | y_{1:T})dx_{t+1}

これにより過去に遡って潜在変数を推定します。

f:id:nogawanogawa:20180721143248j:plain:w500

一括解法

ここからは実際に潜在変数を求めてみます。 まずは過去から知りたい時点までの状態を一括で求める一括解法についてです。

定常な状態時系列(平均・分散・自己共分散・自己相関係数が時間によって変化しない)に対して適用できるウィナーフィルタを見ていきます。

ウィナーフィルタ

下のような観測値と推定値の関係を考えます。

f:id:nogawanogawa:20180915100207j:plain:w500

そして、観測値から状態を推定する逆変換を下のように考えます。

f:id:nogawanogawa:20180915100230j:plain:w500

ここで、


y_t = x_t + v_t


d _t= h_t \bigotimes v_t

となっています。(※\bigotimesは畳み込みを表す)

さて、モデルがなんとなくわかったところで、まずは答えとなる状態と観測値を作ります。


x_t = \phi x_{t-1} + w_t


y_t = x_{t} + v_t

状態は一つ前の状態とノイズw、観測値は状態とノイズvに依存していますね。

この逆変換(ウィナー平滑化)はこんな感じになります。


d_t = \frac{(\phi^{-1} - \beta)(\phi - \beta)}{1-\beta^2} \sum^{\infty}_{k=-\infty} \beta^{|k|}y_{t+k}


\beta = \frac{( \frac{1}{r\phi} + \phi^{-1} + \phi) - \sqrt{( \frac{1}{r\phi} + \phi^{-1} + \phi)^2 -4}}{2}


\gamma = V/W



Σ(゚Д゚)



もはや数学っすね。詳しくはよくわかんないですが、スペクトル変換をして周波数成分に一回落としてからもう一度足し込むみたいです。←適当

これをRで書くとこんな感じになるらしいです。

gist8e32abba01868303601d6ba7cd17c89f

スッキリしてますねー。

こんな感じで観測値から状態を一括推定できました。

感想

難しいところを思い切りすっ飛ばしたのですが、それでも難しいです。 カルマンフィルタに行き着く前に力尽きました。。。

【R】時系列分析の覚書(導入)

こんな本を買ってみたので、エンジニアというよりデータサイエンティストっぽい機械学習の勉強もぼちぼち再開します。

この本すごい評判が良いそうで、読んでて構成がわかりやすかったです。 全体構成と今回の対象は下図にします。

f:id:nogawanogawa:20180715133805j:plain:w500

ということで、まずは導入をやってみます。

時系列分析とは

時系列分析とは、時間方向に関して得られる系列的なデータに対する分析のことです。

もうちょっと詳しく言うとこんな感じらしいです。

基本的には関心のある事象における過去・現在・未来の値を適切に把握し(推定し)、関連してその結果を元に、事象の仕組みや影響に関する知見を得たり対策を考えたりする営みである

やることのなんとなくのイメージがこんな感じです。

f:id:nogawanogawa:20180716101546j:plain:w500

時系列分析を通じて、何らかの値を推定することになるのですが、この本では確定的手法と確率的手法の2種類に分けて解説されています。 確定的手法は状態空間モデルを用いた分析(基本)、統計的手法は状態空間モデルを用いた分析(応用)で扱うようです。

確率・統計について

この本とんでもなくわかりやすいので、詳しい説明は本を読んでいただければと。 用語の定義とかをさわりだけ列挙していきます。

平均やら分散やらは、過去にやっているのでこちらをご参照。

tsunotsuno.hatenablog.com

その他、新しく出てきたとこだけ見ていきます。

複数の確率変数の関係

不安定に揺らぐ変数を確率変数と呼び、大文字のアルファベットで表します。 確率変数自体は値ではなく関数に近い意味合いを持ち、確率変数の具体的な値は実現値と言い、小文字のアルファベットで表します。

確率変数X, Yについて、X, Yが同時に成立する同時確率


p(x, y)

と表されます。 また、確率変数X, Yについて、Yが確定したときのX条件付き確率


p(x | y)

と表されます。 条件付き確率に対して、確率の乗法定理を使用すると以下の関係が導けます。


p(x | y)p(y) = p(x, y) = p(y, x) = p(y | x)p(x)

上式から、


p(x | y) = \frac{p(y | x)p(x)}{p(y)}

が得られます。上の関係をベイズの定理と呼びます。 ちなみに、p(x)事前確率p(x | y)事後確率p(y | x)を(xの条件での)尤度(ゆうど)と呼んだりします。事前確率を条件を確定することで事後確率に変換することをベイズ更新と言ったりします。

確率過程

確率変数が時間方向につながった集合を確立過程といいます。 表記としてはX_t, Y_tのように下付き文字を使用してタイムステップを表現します。

共分散・相関

ある時系列を説明する際に使用される統計量として、共分散・相関が使われます。

共分散Cov)は、確率変数のX, Yの関連を表していて、それぞれの確率変数の期待値と比較したときの値の大小の関係性を表現します。


Cov[X, Y] = E[(X - E[X])(Y-E[Y])]

共分散を規格化したものを相関係数\rho)といいます。


\rho = \frac{Cov[X, Y]}{\sqrt{Var[X]}\sqrt{Var[Y]}}

相関係数は-1〜1の値になり、負、0、正に応じて下記のような関係性があることを示しています。

f:id:nogawanogawa:20180716151100j:plain:w500

確率過程X_t, Y_tに関して、時間軸方向にラグkだけずれたX_t, Y_{t-k}に関する共分散を表す相関関数は、下記のように表現されます。


R(t, k) = E[X_t Y_{t-k}]

異なる確率変数間だけでなく、一つの確率変数の時間方向のラグに関する相関を見る場合は、自己共分散・自己相関係数を使用します。

自己共分散は下記の式で表されます。


Cov[X_t, X_{t-k}] = E[(X_t - E[X_t]) (X_{t-k} - E[X_{t-k}])]

自己相関係数は下記の式で表されます。


\rho_{t,k} = \frac{Cov[X_t, X_{t-k}]}{\sqrt{Var[X_t]}\sqrt{Var[X_{t-k}]}}

自己相関係数も-1〜1の範囲で値を取りますが、今回は値ではなくtkの関数になっています。 時刻とラグの関数から周期的な変動を確認することができます。

定常過程と非定常過程

強定常・弱定常・非定常に関して、時刻が変動したときに平均・分散・自己共分散・自己相関係数・確率分布そのものが変動するかどうかをまとめると下記の様になるそうです。

分類 平均・分散・自己共分散・自己相関係数(○:変化する、☓:変化しない) 確率分布(○:変化する、☓:変化しない)
強定常
弱定常
非定常

この本では、弱定常のことを定常と呼んでいるそうです。 要するに、周期性があれば定常、なければ非定常といったところでしょう。

最尤推定ベイズ推定

確率分布を特徴づけるパラメータ\thetaは一般に未知であり、何らかの方法で特定化します。 確率過程Y_t (t=1, 2, ..., T)全体の尤度はp(y_1, y_2, ..., y_T; \theta)で表されます。 この自然対数である対数尤度


l(\theta) = log p(y_1, y_2, ..., y_T; \theta)

で表され、これを用いて尤度が最大にするように\thetaを決定する方法を最尤法と呼びます。 最尤推定


\hat{\theta}= argmax\ l(\theta)

で表されます。

パラメータ\thetaすら確率変数として取り扱うのがベイズ推定と言うそうです。

とりあえずやってみる

数式ばっかりなのは苦手なので、実際にやってみて感じを掴んでみます。 基本的な分析の流れはこんな感じらしいです。

  1. 目的の確認とデータの収集
  2. データの下調べ
  3. モデルの定義
  4. パラメータ値の特定
  5. フィルタリング・予測・平滑化の実行
  6. 結果の確認と吟味
  7. 1.へ戻る

1. 目的の確認とデータの収集

データ

今回使用するデータはこちらです。

  1. ナイル川の年間流量
  2. 待機中の二酸化炭素濃度
  3. 英国の四半期ごとのガス消費量

※本当は4つ目もあったんですが11章に詳しい説明があるらしく、今回は省略。

目的

分析の目的と対象となるデータを以下に示します。

  目的 データ
1 データ取得中に各時点でのノイズをできる限り除去する 1. 2. 3.
2 過去の急激な変化を捉える 1.
3 未来の値を予測する 2.

2. データの下調べ

まずはデータをそのままの形で見てみます。 データの眺め方にもいろいろあるみたいで、この本では下のような表示を行っています。

横軸時間のプロット

ナイル川の流量については、あまり規則性が見られません。 これは不規則な擾乱を除けば、毎年概ね同じ値に保たれているということを意味します。 二酸化炭素濃度については、小刻みに上下を繰り返しながら、全体としては右肩上がりの傾向が見て取れます。 イギリスのガスの量は右肩上がりの傾向があるようで、更に年々変動の幅が大きくなっていることがわかります。

ヒストグラムと五数要約

どのグラフにも大きな外れ値はなさそうです。 ここで、二酸化炭素濃度のデータについては、東日本大震災時の際のデータの欠損がありますので、前後の月のデータの値の平均で代用しています。 欠損は放って置くと変な結果になったりするそうなので、潰しておきます。

自己相関係数

全体として自己相関係数の値がある程度ありそうなので、時系列的な関係性がありそうだとわかります。 また、イギリスのガスのデータはラグが整数のとき自己相関係数が大きくなる傾向があるようです。

周波数スペクトル

フーリエ変換して周波数成分の大きさを確認します。 音とかで言うところの、「聞こえてきた音の中から、成分("ド"とか"レ"とか)がどれくらい入っているか」を調べる作業みたいなものです。

ナイル川の流量についてはあまり周波数成分に特徴が見られません。 二酸化炭素濃度とイギリスのガスについては、特定の周波数成分が大きくなっている事がわかります。

3. モデルの定義

だいたい、傾向がわかったところで仮置きでモデルを設定してみます。

f:id:nogawanogawa:20180717213549j:plain

左から、レベル、傾き、周期を表していて、これを組み合わせでグラフが構成されていそうだということが想像つきます。

  1. レベル
  2. レベル + 傾き + 周期
  3. レベル + 傾き + 周期

なんとなく、こんな感じの組み合わせで表現できそうです。

4. パラメータ値の特定

検討したモデルに対して、具体的な数字を当てはめていきます。 上で指定したモデルでは、経済分野などで使用されるホルト・ウィンタース法が使用できます。

ホルト・ウィンタース法

k期予測値\hat{y}_{t+k}は次のように表されます。


\hat{y}_{t+k}= level_t + trend_t k + season_z{t-p+k_p^+}

このとき、レベル(level)、傾き(trend)、周期(season)はそれぞれ、

level_t = \alpha(y_t - season_{t-1}) + (1 + \alpha)(level_{t-1} + trend_{t-1})

 trend_t = \beta(level_t + level_{t-1}) + (1-\beta)trend_{t-1}

season_t = \gamma(y_t - level_t) + (1 - \gamma)season_{t-p}

となります。ここで、k_p^+ =  \lfloor(k-1)mod p\rfloor +1です。

式で書くとややこしいですが、Rにはこれらをビシッと計算してくれる関数が用意されており、HoltWinters()でパラメータを決定できます。

5. フィルタリング・予測・平滑化の実行

HoltWinters()でサクッとパラメータを特定できます。

赤い点線がモデルによるフィッティング(フィルタリング値)を表すグラフです。

成分ごとに見ていくとこんな感じだそうです。

未来の予測に関しては、二酸化炭素濃度を使用して、未使用の2015年分を予測するとこんな感じになります。

6. 結果の確認と吟味

最後に仮定したモデルが良いか悪いかを評価します。 今回、目的はこんな感じでした。

  • データ取得中に各時点でのノイズをできる限り除去する
  • 過去の急激な変化を捉える
  • 未来の値を予測する

ノイズの除去に関しては、残差(Residual)を確認します。これが大きいとモデルと実測値がかけ離れていることになるので、ノイズを除去しきれていないことを意味します。これはどこまでを許容するかわかりませんが、今回は一旦良しとします。

また、過去の急激な変な変化を捉えるという点については、ナイル川の流量が1899年に急減しているそうですが、そこを捉えることはできていません。この点については更に改良が必要と言うことがわかります。

未来の予測については、二酸化炭素濃度の予測は直感的には一致しているように見えるので、問題ないこととします。

ここで出てきた問題を踏まえて再度手順をやり直し、モデルを改良していきます。

感想

基本的な用語、式の説明と簡単な時系列分析のやり方の説明でした。 久しぶりにやるとRむずいっす。。。

【論文メモ:HDGAN】Photographic Text-to-Image Synthesis with a Hierarchically-nested Adversarial Network

論文

[1802.09178] Photographic Text-to-Image Synthesis with a Hierarchically-nested Adversarial Network

著者

Zizhao Zhang, Yuanpu Xie, Lin Yang
University of Florida

背景

テキストから写実的画像を生成すること(Text-to-Image)は生成モデルの分野での大きなテーマとなっている。 このテーマでは、精細な画像を生成することと入力したテキストとの一貫性を両立することが必要になる。

Text-to-Imageの分野では、近年ではGANを用いた生成モデルが多数提案されている。 GANを使用する場合には学習の収束と安定性が大きな課題となる。 これらの課題を解決し解像度とテキストの一貫性を実現するような先行研究が行われてきたが、解像度が低かったり大きなネットワークを必要としたりと、まだまだ確立した手法があるとは言えない。

目的とアプローチ

目的

高解像度かつテキスト一貫性を併せ持つ生成モデルの実現

アプローチ

HDGAN (High-Definition results and the idea of Hierarchically-nested Discriminators)

  • one generator - multi discriminator の使用によるステージングを必要としない
  • 中間生成画像を段階的にDiscriminatorに入力するLoss

提案手法

HDGAN

GANを使用した生成モデルの先行研究では、学習を安定させるためにStageを使用したり、multi Discriminatorにすることで、DiscriminatorがGeneratorより過剰に学習していくことを回避していた。 今回の提案手法であるHDGANの概要を下記に示す。

f:id:nogawanogawa:20180625223238j:plain

HDGANは、Generator1つに対して、複数の階層的Discriminatorを用意する。 Generatorの中間生成画像をそれぞれDiscriminatorに入力することで、Generatorを1つにすることが可能である。

また、先行研究と提案手法を比較を下記に示す。

f:id:nogawanogawa:20180625223635j:plain

(A)StackGANでは、低解像度の画像を生成するStage-Ⅰと高解像度の画像を生成するStage-Ⅱの二段階のステージを使用する。 Stage-Iを先に重点的に学習することで、安定した学習が可能になる一方で、通常のネットワークに加えステージを導入するためアーキテクチャが複雑になる。

(B)また、GMANのようにDiscriminatorを複数使用したGANのモデルも提案されている。 これによって、Discriminatorに1つの場合より、Generatorの学習が安定する。 一方で、これを高解像度の場合に適用してうまくいくかは不透明である。

(C)さらに、PGGANでは低解像度から高解像度まで段階的にGeneratorとDiscriminatorの層を増やしていく。 これにより、低解像度の層が安定するまで学習して層を追加するので、安定して学習させる事ができる。 一方で、高解像度を生成しようとするにつれて、高い演算性能ができる。

(D)提案手法では、Generatorの中間生成物と対応する解像度のDiscriminatorに入力する。 これにより、隠れ層の出力に対するフィードバックの信号がより強いフィードバックとして反映されるため、学習が安定化する。

Multi-purpose adversarial losses

一般的なGANの目的関数を下記に示す。

f:id:nogawanogawa:20180702213014j:plain

ここでGはGenerator、DはDiscriminatorを表している。 このように一般的なGANではGeneratorとDiscriminatorのLossを最小化することが目的関数になっている。 このときのLossは下記の様になっている。

f:id:nogawanogawa:20180702213741j:plain

Eはクロスエントロピーを表し、DiscriminatorをGeneratorに対して2つの項を考えている。

一方、HDGANの目的関数を下記に示す。

f:id:nogawanogawa:20180702213105j:plain

一般的なGANに対して、Textに関する項tが追加されている。 このときのLoss計算の概念図とLossの計算式を下記に示す。

f:id:nogawanogawa:20180625223709j:plain

f:id:nogawanogawa:20180702213118j:plain

HDGANでは、Discriminatorが複数あるmulti-Discriminatorの構成を取るため、Generatorの隠れ層の分だけLossを個別に計算して足し合わせ る。 そのため、Lossも各階層について計算した後、総和を取ることで表現される。

Architecture Design

Generator

Generatorは単純な3つのモジュールで構成される。

  • K-repeat Res block
    • 畳み込み層×2
    • ReLU
    • batch normalization
  • stretching layers
    • サイズ2の近傍アップサンプリング
    • ReLU
    • 畳み込み層×1
    • batch normalization
  • linear compression layers
    • 畳み込み層
    • Tanh(RGB空間に圧縮)

入力は conditioning augmentationと学習済みのembedding matrixによって生成された1024×4×4のembeddingを使用する。

Discriminator

Discriminatorはシンプルにストライド2の畳み込み層に対してバッチノーマライゼーションとLeakyLUを適用させる。 実装には2通りの選択肢がある。 1つはDiscriminatorそれぞれに対して教師データと生成データの判別を行う方法である。 もう1つは、始めに512×4×4の特徴マップと128×4×4のtext embeddingを合体させる。 それからサイズ1の畳み込み層を使用してテキストと画像の特徴を統合し、4×4の畳み込み層によって画像とテキストのペアに対して教師データと生成データとを判別する。

その他の最適化等は通常のGANと同様であるため割愛する。

評価

Experimental Setup

Dataset

評価に使用するデータセットは下記の3種類である。

  • CUB dataset
  • Oxford-102
  • COCO dataset

Evaluation metric

評価の指標としてInception scoreとMulti-scale structural similarity (MS-SSIM score)を使用する。 これら2つはGANの評価に広く使用されている一方、テキストと生成画像の一貫性については評価できない。 そこで本研究では、 Visual-semantic similarityを評価指標として導入する。

Visual-semantic similarityは下記の式によって評価する。

f:id:nogawanogawa:20180702213933j:plain

ここでvはインセプションモデルによって抽出された特徴ベクトルを表す。 スコアリング関数cを使用して二次元のロスを計算している。

テスト段階においては、text embeddingと生成画像のペアで算出される特徴ベクトルは同じであることが望ましい。 上の式によって、ペアの特徴ベクトルが親しいほど、大きく算出される仕組みとなっている。

先行研究との比較

提案手法の妥当性を評価するために、先行研究との比較を行う。 比較対象は、GAN-INT-CLS、GAWWN、TAC-GAN、Progressive GAN、StackGAN、StackGAN++をとする。 特に、StackGANについては詳細に評価する。

提案手法と比較対象のInception scoreを下記に示す。

f:id:nogawanogawa:20180625224100j:plain

HDGANでは、CUBのデータセットでStackGANに比べて45%、StackGAN++に比べて31%良いなど、先行研究に対して非常に良い効果が得られている事がわかる。 Oxford-102では、TAC-GANと同等の結果が得られているが、TAC-GANで使用しているような付加情報を本研究では使用していないという優位性がある。

提案手法とStackGANで生成された画像を下記に示す。

f:id:nogawanogawa:20180625224207j:plain

f:id:nogawanogawa:20180625224241j:plain

提案手法のほうが、文章の詳細も表現され、自然な色使いで複雑な形状を表現できている。

解像度別のInception scoreを下記に示す。

f:id:nogawanogawa:20180625224313j:plain

解像度が高くなってもInception scoreはStackGANより優れている。 また、StackGANでは画像の一貫性が崩れているが、HDGANでは解像度を大きくしたときに画像の一貫性が崩れていない。

Visual-semantic similarityの評価を下記に示す。

f:id:nogawanogawa:20180625224124j:plain

HDGANのほうが値が大きく、意味論的な一貫性が取れていることがわかる。

さらに多様性についても評価した結果を下記に示す。

f:id:nogawanogawa:20180625224438j:plain

HDGANは、同じ入力からでも広いバリエーションで画像を生成できている。

MS-SSIMについてStackGAN・Prog.GANと比較した結果を下記に示す。

f:id:nogawanogawa:20180625224507j:plain

StackGANよりEquallity lineより小さく、多様性の高い優れたモデルとわかる。 また、Prog.GANと比較してもMS-SSIMの値は優れている。

文を書き換えによるstyle transfer

文の書き換えによるstyle transferに関する結果を下記に示す。

f:id:nogawanogawa:20180625224528j:plain

提案手法により、なめらかに画像が変換されている事がわかり、細部まで表現されている。

考察

階層的敵対性学習

階層化した学習の妥当性に関して、下記に示す。

f:id:nogawanogawa:20180625224558j:plain

Discriminatorの層が深くなればなるほど、当手法は有効であると考えられる。 また、StackGANではStageごとにTextを入力していたが、当手法では入力は一度だけにもかかわらず、良い結果となっている。

局所的画像のLoss

局所的画像のlossの使用について評価する。

f:id:nogawanogawa:20180625224614j:plain

上の表ではlocal image lossを使用した場合(w/)、使用しない場合(w/o)よりも値が優れることがわかる。 そのため、local image lossの活用は妥当であるといえる。

また、下記の図からも、テキストの内容をより詳細に反映した画像を生成できることがわかる。

f:id:nogawanogawa:20180625224632j:plain

結論

本論文では、中間生成画像をDiscriminatorに階層的に学習させるHDGANを提案した。 本手法により、先行研究よりテキストに忠実かつ表現力のある生成モデルが実装されることがわかった。

【論文メモ:PGGAN】Progressive Growing of GANs for Improved Quality, Stability, and Variation

論文

https://arxiv.org/abs/1710.10196

著者

Tero Karras, Timo Aila, Samuli Laine, Jaakko Lehtinen
NVIDIA

背景

生成的手法の中でも現在特に優れているものには、autoregressive models、variational autoencoders (VAE) そして generative adversarial networks (GAN)などが挙げられる。 これらにはそれぞれ長所・短所があり、現在も研究が進められている。

このうち、GANは生成する分布(画像など)の解像度が高くなるに連れてランダム要素の影響が色濃くなり、DiscriminatorはGeneratorの生成分布と教師データを区別することが容易になってしまい、ネットワーク全体の学習が不安定になる。 このように、GANは高解像度の画像を出力することが難しいという問題がある。

目的とアプローチ

目的

  • GANを用いた高解像度分布の生成

アプローチ

PGGAN(Progressive Growing of GANs)

  • 段階的にネットワーク層を増加させ、追加した層の影響を学習に度合に応じて比重を変化
  • ミニバッチ標準偏差を活用した多様性の向上
  • マルチスケールによるワッサースタイン計量を用いた統計的類似度の使用

提案手法

PGGAN(Progressive Growing of GANs)

PGGAN(Progressive Growing of GANs)の概念図を下記に示す。

f:id:nogawanogawa:20180528215414j:plain

PGGANでは、低解像度の画像(上では4×4)から学習を開始し、徐々に層を追加することで対象の画像の解像度を向上させる。 層を新しく追加する際には、学習がなされていないためランダム要素が強く、Generatorの出力にランダム要素が強く反映されてしまう。 そこで、PGGANでは下図のように、出力に対する新しく追加した層の影響度(α)を段階的に増やす。

f:id:nogawanogawa:20180528215600j:plain

これにより、ネットワーク全体の出力として学習が完了していない層のランダム要素を軽減することが可能になり、学習が安定化する。

ミニバッチ標準偏差を活用した多様性の向上

GANでは、学習データの範囲でしか特徴を学習できず、学習データの限られた特徴に大きく影響された学習をしてしまう恐れがある。 この問題を解決するために、ミニバッチの標準偏差を用いた正規化によってミニバッチのまとまりから得られる特徴を学習する方法を導入する。

はじめに、ミニバッチ内のテンソル(H×W×C)で同じ位置に対応する画素について標準偏差を計算し、H×W×Cのテンソルを一つ求める。 次にそのテンソル内のすべてのピクセルの平均を計算し、スカラ値を一つ得る。 最後にそのスカラ値を複製し、一つのテンソルを作成する。

この処理は、理論上ネットワークのどこにでも挿入することができるが、諸々の事情からDiscriminatorの最終層に挿入する。 (詳細は論文のAppendix参照)

補足

イメージとしてはこんな感じらしいです。

f:id:nogawanogawa:20180603112433j:plain

GeneratorとDiscriminatorの正規化

学習データに特徴的な学習データが含まれている際に、そのデータに学習結果が強く影響を受けてしまう問題がある。 一般にこの問題に対してはバッチノーマライゼーションを使用して対処することが多い。 しかし、我々はGANのネットワーク自体が問題なのではなく、学習信号の大きさとその評価に制限がないことに問題があると考えた。

そこで我々は下記の2つの成分を使用することでこの問題に対処する。

Equalized learning rate

従来のネットワークの重みの繊細な初期化は行わず、0~1の範囲での初期化に対して実行時にネットワークの重みを変化させる。 ネットワークの重みに下記の式を導入する。


\hat{\omega_i} = \omega_i / c

ここで、\omega_iはネットワークの重み、cはレイヤごとの標準化定数を表す。

通常Adamなどでは、入力パラメータとは独立に標準偏差をもとに勾配を決める。 そのため、ダイナミックレンジが大きいパラメータはそうでないパラメータよりも学習に時間を要する。 我々の手法では、ダイナミックレンジを考慮されるため、すべてのパラメータで均質の学習スピードを確保することができる。

Pixelwise feature vector normalization in generator

GeneratorとDiscriminatorがコントロールできなくならないように、Generatorの各畳み込み層のあとに各ピクセルごとにfeature vectorを正規化する処理を行う。

f:id:nogawanogawa:20180603122733j:plain:w280

ここで、\epsilon=10^{-8}はネットワークの重み、Nは特徴マップの数(チャンネル数)、abはそれぞれオリジナルと正規化された特徴ベクトルを表す。

これにより、変化の小さな学習の際には影響は小さく、変化の大きな学習の際には必要に応じて信号を減衰させることができる。

マルチスケールにおける統計的類似性の使用

異なるGANのモデルを評価するには、膨大な数の生成画像を評価する必要があり、主観的になってしまう。 したがって、これらの評価には画像から得られる指標を用いて自動的に評価する必要がある。 これまでも、multi-scale structural similarity (MS- SSIM)のような手法もあったが、色や形状の多様性による細かな効果には影響しにくい。

我々は、優れたGeneratorはどんなスケールにおいても学習データと類似した局所的画像構造を生成するものだと考えている。 そこで、教師データとラプラシアンピラミッドによって得られる局所イメージを統計的に比較する手法を用いる。 連続するピラミッド階層において、前のレベルをアップサンプリングしたものと次のレベルとの差分を見ることで、これを実現する。

各層は7×7の画素近傍によって表現される記述子(各層128個)で表現される。 それぞれの記述子を正規化し、標準偏差を計算、その後ワッサースタイン計量によって類似度を判定する。 これにより、小さい解像度部分で見たときの差分は画像全体の大まかな構造を捉え、高解像度部分で見たときの差分は輪郭やノイズといった部分まで表現する。

補足

イメージとしてはこんな感じみたいです。

f:id:nogawanogawa:20180603153152j:plain

評価

実際に生成された画像はこちら。

youtu.be

統計的類似度に関する評価

sliced Wasserstein distance (SWD) とmulti-scale structural similarity (MS- SSIM)に関する評価を行う。

条件は下記のとおり。

  • 最新のloss function (WGAN-GP)と学習の設定を使用
  • 学習データ : CelebA, LSUN BEDROOM
  • 解像度 : 128×128

次の表は数値的にSWDとMS-SSIMに関してワッサースタイン計量とMS-SSIMを計算した結果である。

f:id:nogawanogawa:20180528215914j:plain

また、提案手法によって提案された画像を下記に示す。

f:id:nogawanogawa:20180528220113j:plain

直感的に、良い評価指標では様々な色や形状、視点が反映されたもっともらしい画像が良い評価とされるべきである。 つまり、(h)の行が最も値が小さくなるはずである。 しかし、今回それがMS-SSIMには見られない。 そのため、MS-SSIMより提案手法のSWDのほうが良い指標と考えられる。

収束及び学習スピード

収束と学習速度に関する評価結果を下記に示す。

f:id:nogawanogawa:20180528220454j:plain

図(a)はprogressive growingがない場合、(b)はprogressive growingがある場合を示している。 これらを比較すると、progressive growingを使用したほうが収束時のSWDの値が小さく、学習時間を短縮できている事がわかる。

progressive growingなしだとマクロな視点とミクロな視点の両方の学習を同時にこなす。 一方progressive growingありの場合は低い解像度の学習は終わっているため、解像度を大きくしたときに安定して学習できている。

図(b)より、各ラプラシアンピラミッドの層において、はじめの低解像度のときはまっすぐSWDが減少している。 解像度が向上しても、一貫してSWDは減少傾向を示している。 一方、図(a)ではすべての各ラプラシアンピラミッドの層で似たカーブの動きを見せ、データによって共鳴した学習をたどっている。

図(c)より、progressive growingを使用した場合は1024×1024の画像について収束するまでに96時間要したのに対し、progressive growingを使用しなかった場合には520時間を要した。 そのため、progressive growingを使用したほうが約5.4倍高速化する事がわかる。

CelebAを用いた高解像度画像の生成

本論文で1024×1024の高解像度の画像を生成した結果を下記に示す。

f:id:nogawanogawa:20180528220626j:plain

これらの画像を生成するのに、Tesla V100 GPUを8機使用して、4日計算した。

LSUNの結果

LSUNのBED-ROOMについて、生成画像を比較したものを下記に示す。 f:id:nogawanogawa:20180528220756j:plain

また、下記には別のLSUNカテゴリの図を示す。

f:id:nogawanogawa:20180528221024j:plain

全般的に生成画像の質は高いが、先行研究と大差ない。

実装

著者の方が実際に書いてくれているので、こちらをご参照。

https://github.com/tkarras/progressive_growing_of_gans

コードを眺めるだけならいいけど、実際に計算は勘弁してください。 論文の実装では100万くらいのGPU8機使ってるし、GitHubのREADMEによるとGPUは1000万超えるみたいですし。 ここまでくると一般人では、予算的に実行することすら不可能ですね。

おまけ

書いてる途中で超わかりやすいブログで見つけました。 すごい参考にさせていただきました。 ありがとうございます。

st-hakky.hatenablog.com

TensorBoardを使ってみた

今回はTensorFlowのライブラリに付属しているTensorBoardなるものを使って色々可視化してみたいと思います。

参考にさせていただいたのはこちらです。

deepage.net

TensorBoard

TensorFlowにはTensorBoardという付属のツールがあって、こいつを使えるのがTensorFlowの結構大きなポイントだったりします。

TensorBoardではTensorFlowでの実行時に発生するデータを可視化・表現できます。 可視化できるものとして代表的なものはこんな感じらしいです。

使い方

TensorBoardでは、可視化したいデータや形式をプログラム中に明示的に指定することで、 TensorBoardが表示するための情報を裏で取得してくれるみたいです。

手順

そもそも動くTensorFlowを使用したコードがある状態から考えると、手順としてはこんな感じでしょうか。

  1. プログラム中に可視化用のコードを追加
  2. 普通にプログラムを実行
  3. TensorBoardを起動
  4. ブラウザからTensorBoardを開く

具体例

今回は試しにこの前使用したコードを使用してやってみます。

tsunotsuno.hatenablog.com

1. プログラム中に可視化用のコードを追加

使うためには、上のコードにちょびっと追加します。

# TensorBoardで追跡する変数を定義                                                                                                                               
with tf.name_scope('summary'):
    tf.summary.scalar('loss', cross_entropy)
    merged = tf.summary.merge_all()
    writer = tf.summary.FileWriter('./logs', sess.graph)

2. 普通にプログラムを実行

特に何も考えずに普通に実行してください。

python3 main.py

普通に実行します。 すると、カレントディレクトリにlogというディレクトリができているかと思います。 これができていれば、とりあえず問題なく動いていそうです。

3. TensorBoardを起動

実行が完了したら、TensorFlowを起動します。 TensorBoardは次のコマンドで起動します。

tensorboard --logdir=./logs

4. ブラウザからTensorBoardを開く

コンソールに出ているURLをブラウザに放り込めばTensorBoardを確認する事ができます。 きちんとTensorFlowが起動できていれば下のurlをで可視化結果が開きます。

http://localhost:6006/

f:id:nogawanogawa:20180131203555p:plain

Google Cloud Platformの導入メモ

この前はGoogle Colabを使って見たんですが、どうやらGANをやるには不向きな感じでした。

tsunotsuno.hatenablog.com

こうなったら、やはり本家だろうということでGoogle Cloud Platform (GCP)を使って機械学習をやる方法を調べてみました。

参考にさせていただいたのはこちら。

qiita.com

Google Cloud Platform (GCP)

登録

何はともあれ登録。

cloud.google.com

料金

無料トライアルで始めの1年は$300分使えるみたいです。 (2018年5月3日現在)

f:id:nogawanogawa:20180503223011j:plain

その他の料金はこちら。

cloud.google.com

クラウド使うのが初めてなので、どれ使えばいいとかイマイチわかりません。 参考サイトでも言われていますが、すぐに$300使い切ることはまず無いので、一旦使ってみてからサービスは考えることにしたほうが良さそうです。

チュートリアル(途中まで)

Google公式の説明書はこちら。

(英語)https://cloud.google.com/ml-engine/docs/tensorflow/getting-started-training-prediction?hl=ja

(日本語)https://cloud.google.com/ml-engine/docs/getting-started-training-prediction

どうやら、CLOUD SHELLを使うほうがお手軽らしいので、今回はそちらでやってみました。

"TensorBoard を使用した要約ログの検査"までやってみましたが、書いてあるとおりにコマンドをコピペしていけばいけました。

f:id:nogawanogawa:20180504100030p:plain

f:id:nogawanogawa:20180504100039p:plain

あとはコードを手元のコードを吸い上げれば普通に実行できそうです。

まとめ

GCPチュートリアルレベルであれば結構簡単に動きました。 クラウドとのファイル転送やらGPUの使い方やらは勉強しないといけませんが、ひとまず環境が使えることを確認できてよかったです。