「階層ベイズモデルとその周辺」を読む～第II部非線形ダイナミカルシステムの再構成と予測～

はじめに

本書は，岩波書店「統計科学のフロンティア」シリーズ(全12巻)の第4巻である。初版発行は2004年であり，当時最先端だった統計科学の話題を，第一人者の先生方がコンパクトにまとめているというものである。

第4巻「階層ベイズモデルとその周辺―時系列・画像・認知への応用―」は3部+補論の構成となっており，以下のような章立てになっている。

第I部　事前情報を利用した複雑な系の解析　(石黒真木夫)
第II部　非線形ダイナミカルシステムの再構成と予測　(松本隆)
第III部　視覚計算とマルコフ確率場　(乾敏郎)
補論　帰納推論と経験ベイズ法―逆問題の処理をめぐって―　(田邉國士)

本記事は，第II部　非線形ダイナミカルシステムの再構成と予測の読書メモである。

www.iwanami.co.jp

はじめに
- 目次
本書・第II部の興味深かった点
第1章　問題提起と導入
- ダイナミカルシステム
- 非決定論的ダイナミカルシステム
第2章ニューラルネットワーク
第3章ダイナミカルシステムの学習と予測
- 3.1 モデル定式化
- 3.3 予測アルゴリズム
第4章具体的問題
第5章ハミルトニアン・モンテカルロによるベイズ的学習と予測
本書に関するまとめと感想
- モデル選択も含めて階層ベイズモデルの枠組みで行なっていることが参考になった
- 実装にチャレンジしたい

本書・第II部の興味深かった点

本書・第II部のテーマは，非線形ダイナミカルシステム(ノイズを含む非線形な時系列モデル)を，ニューラルネットワークでモデリングし，階層ベイズモデルの枠組みでデータに当てはめる，というものである。

本書において，興味深かった点は以下の3点である。

1. モデル選択も含めて階層ベイズモデルの枠組みで行なっている。

時系列モデルを，ニューラルネットワークを用いて推定するので，時間遅れ $\tau$ や，中間層の素子数 $h$ といったハイパーパラメータが出てくる。
このハイパーパラメータについて，グリッドサーチを用いることなく，「モデル周辺尤度」という枠組みを用いてベイズモデルの枠組みで完結させている。

2. 事後分布を様々な近似を用いて求めている。

ベイズ統計で難しいのは，事後分布や予測分布を求めることである。本書でも予測分布を求めることがテーマの1つであるが，予測分布を2次近似やMAP推定量を駆使することで近似的に求める方法を説明している。
これにより，予測分布の平均値を求める方法を説明している。

3. ハミルトニアン・モンテカルロ法を用いた事後分布推定の方法を説明している。

ベイズ統計学の強みは，予測値を点で求めるのではなく，幅を持った分布で求めることができる点である。事後分布を求める方法として，マルコフ連鎖モンテカルロ法(MCMC)の1種であるハミルトニアン・モンテカルロ法を用いる方法を紹介している。

第1章　問題提起と導入

本章では，分析対象とする非線形ダイナミカルシステムについて説明している。

ダイナミカルシステム

ダイナミカルシステムとは，ある時点 $t+1$ の値 $x_{t+1}$ が，過去の値 $x_t,...$ によって決まるシステムのことである。

$\begin{align} \text{離散系 : } x_{t+1} &= F(x_t), \quad t=0, 1, 2, \cdots \\ \text{連続系 : } \frac{d}{dx} x_{t} &= F(x_t) \\ \end{align}$

本書では，このうち離散系について説明する。

関数 $F$ が線形であれば，将来値である $x_{t+1}$ の挙動(解軌道)を予測することは扱いやすいが，非線形の場合，たとえば

$\begin{align} x_{t+1} = a x_t(1 - x_t) \end{align}$

のようなものでも解軌道は複雑になることが知られている。

また，上記のように入力を含まない系を自律ダイナミカルシステム，入力 $u_t$ を含むシステム $x_{t+1} = F(x_t, u_t)$ を非自律ダイナミカルシステムと呼ぶ。

非決定論的ダイナミカルシステム

ダイナミカルシステムのうち，ノイズを含むものを非決定論的ダイナミカルシステムと呼ぶ。
線形のシステムのうち，代表的なものは線形状態空間モデルである。

$\begin{align} x_{t+1} &= Fx_t + G\nu_t \\ y_t &= Gx_t \\ \end{align}$

線形でない非決定論的ダイナミカルシステムは非線形な関数 $f$ を用いて以下のように表される。

$\begin{align} x_{t+1} = f(x_t, x_{t-1}, \cdots, x_{t-\tau+1}; w) + \nu_t \\ \tag{1} \end{align}$

非線形性は，非線形な基底関数の組合せによって表現する。また $w$ はデータから求められるパラメータ(回帰分析における回帰係数に相当)である。
本書の分析対象は，(1)式で表されるような式である。

第2章ニューラルネットワーク

本書では，非線形な関数として，ニューラルネットワークを用いる。ニューラルネットワークモデルを用いる際に検討することは以下の3つである。

関数族を選ぶ
1. ニューラルネットワークは，シグモイド関数や動径基底関数といった非線形な基底関数の組合せからなる。分析者は，この関数族を選ぶ。
パラメータを調整する
1. データにフィットするように，パラメータの学習を行なう。
パラメータを選ぶ
1. これは，ニューラルネットワークにおける中間層の素子数などに相当するが，これを選択する。

第3章ダイナミカルシステムの学習と予測

本節では，ダイナミカルシステムの予測を，ベイズ統計の枠組みとしてとらえて，ニューラルネットワークの学習(再構成)と予測を行なうための方法を説明している。
なお，非線形ダイナミカルにおける不確定性(システム雑音，初期値，観測雑音)のうち，システム雑音のみを対象とする(すなわち，初期値は観測可能で，状態量もそのまま観測できる，とする)。

3.1 モデル定式化

非線形ダイナミカルシステムを対象に行ないたいことは，過去の時系列データ $\{ x_t\}, t=0, 1, \cdots, N$ が与えられたときに，将来値 $\{ x_t\}, t \gt N$ を予測することである。

非線形ダイナミカルシステムを統計モデルに落とし込むと，下図のようになる。

この統計モデルの各要素の詳細を説明する。

(i)アーキテクチャ

ユーザが決めておくモデルの構成要素であり，

ニューラルネットの基底関数 $f$ 　(シグモイド関数，動径基底関数など)
時系列モデルの時間遅れ $\tau$
ニューラルネットワークの隠れ層の素子数 $h$

などである。
なお機械学習の枠組みでは，これらの最適化はグリッドサーチで行なわれるが，本書ではモデル周辺尤度という概念を導入し，ベイジアンモデリングの枠組みで最適化を行なう。

(ii)尤度

非線形ダイナミカルシステムの予測モデルを

$\begin{align} x_{t+\tau} = f(x_{t+\tau-1}, \cdots, x_t; w) + \nu_{t+\tau}, \quad \nu_{t+\tau} \sim \mathcal{N}(0, 1/\beta) \\ \end{align}$

とする。

学習データ $D = \{ x_0, \cdots, x_N \}$ について尤度は

$\begin{align} P(D| w, \beta, H) &= \prod_{t=0}^{N-\tau} \frac{1}{Z_D(\beta)} \exp \left( -\frac{\beta}{2} (x_{t+\tau} - f(x_{t+\tau-1}, \cdots, x_t; w) )^2 \right) \\ &\times P(x_{\tau-1}, \cdots, x_0 | H) \\ \end{align}$

と書ける。ただし， $Z_D(\beta)$ は規格化定数である。

(iii)事前分布など

階層ベイズモデルを扱うので，パラメータ $w$ にも事前分布を設ける。パラメータ $w = (w_1, \cdots, w_k)^T$ を $C$ 個の部分ベクトルに分割し，

$\begin{align} w = (w_1, \cdots, w_C), \quad w_c \in R^{k_c} \\ \end{align}$

とし，各部分ベクトルパラメータ $w_c$ が正規分布に従うと仮定する。

$\begin{align} P(w | \alpha, H) &= \prod_{c=1}^{C} \frac{1}{Z_w(\alpha_c)} \exp \left( -\frac{\alpha_c}{2} \lVert w_c \rVert ^2 \right) \\ \alpha &= (\alpha_1, \cdots, \alpha_C), \quad \alpha_c \in R \\ \end{align}$

通常のベイズ線形回帰では，パラメータ $w$ の事前分布に設けるパラメータは単一のものとしている(すなわち $\alpha_1 = \cdots = \alpha_c = \alpha$ )。
このように， $w_c$ の値に応じて事前分布のパラメータ $\alpha_c$ を変えることによって，予測に寄与しない $w_c$ に対応する $\alpha_c$ は値が大きくなり，正規分布の分散が小さくなって，値が0に集中する。
すなわちスパース回帰と同様に，不要なパラメータが除去される効果が得られる(関連度自動決定(Automatic Relevance Determination; ARD))。

また， $\alpha$ 自体の事前分布や，尤度に出てくる精度 $\beta$ については無情報なので，一様分布またはガンマ分布を考える。

アーキテクチャに関する事前分布 $P(H)$ にも同様に一様分布を設定する。アーキテクチャに関する事前分布は，例えば $\tau = 1, 2, ...$ のように $\tau$ の値を何パターンか用意しておき，それらについてデータの当てはめを行なうことを考えると，これらのモデルの選び方には差をつけずに選択をするので，一様分布になると考えられる。

グラフィカルモデルによる表現

本書には出てこないが，グラフィカルモデルを用いて変数間の関係を表現してみる。上記の統計モデルから，変数間の関係を表現すると以下のようになる。

(i)アーキテクチャ
- アーキテクチャには，時間遅れの変数 $\tau$ などが含まれる。
(ii)尤度
- 予測モデルでは，データ $D={x_t}$ を，アーキテクチャ $H$ に関するパラメータである $f, \tau$ とパラメータ $w$ ，また精度 $\beta$ で表現している。
(iii)事前分布など
- パラメータ $w$ は $\alpha$ を含む事前分布で表現される。

そのためグラフィカルモデルは以下のように表現される。

3.3 予測アルゴリズム

パラメータの事後分布

上記のモデルにおいて，推定するべきパラメータは以下の5つである。

連続パラメータ
1. ニューラルネットワークのパラメータ $w$
2. 事前分布のパラメータ $\alpha$
3. システムノイズの精度 $\beta$
離散パラメータ (まとめてと表現する)
1. 時系列モデルの時間遅れ(マルコフ過程の次数) $\tau$
2. ニューラルネットワークの中間素子数 $h$

事後分布は，ベイズ公式から以下のように得られる(グラフィカルモデルも参照のこと).

$\begin{align} P(w | D, \alpha, \beta, H) &= \frac{P(D|w, \beta, H) P(w| \alpha, H) }{ P(D|\alpha, \beta, H) } \\ P(\alpha, \beta | D, H) &= \frac{P(D | \alpha, \beta, H) P(\alpha, \beta | H)}{ P(D|H) } \\ P(H | D) &= \frac{P(D | H) P(H)}{ P(D) } \\ \end{align}$

事後確率最大化によるパラメータ推定

連続パラメータ $w, \alpha, \beta$ は，事後確率最大化によって求める。このように求められたMAP推定量を $w_{MP}, \alpha_{MP}, \beta_{MP}$ と表現する。

モデル周辺尤度によるハイパーパラメータの選択

離散パラメータ $\tau, h$ などをまとめて $H$ と表現する。 $H$ も事後確率最大化によって求める，というのが本書の特徴の1つである。

$\begin{align} P(H | D) &= \frac{P(D | H) P(H)}{ P(D) } \propto P(D | H) P(H) \\ \end{align}$

であるので，

$\begin{align} \mathrm{argmax} _ H P(H | D) = \mathrm{argmax} _ H P(D | H) P(H) \end{align}$

となる。
ここで，

$\begin{align} P(D | H) = \int P(D | \alpha, \beta, H) P(\alpha, \beta | H) d\alpha d\beta \\ \end{align}$

であり，これはモデルに関する周辺尤度とみなせる。

実際は，この積分計算をするのが大変である。そこで， $\alpha, \beta$ の分布がそれぞれ $\alpha_{MP}, \beta_{MP}$ に集中していると仮定すると，もはや $\alpha, \beta$ は分布ではなく定数 $\alpha_{MP}, \beta_{MP}$ とみなせるようになり，積分が外れる。また $P(H)$ は一様分布であると仮定しているので，

$\begin{align} \mathrm{argmax} _ H P(H | D) = \mathrm{argmax} _ H P(D | \alpha_{MP}, \beta_{MP}, H) \end{align}$

のように近似することができる。

予測分布

新たなデータ $\{ x_{N+t} \}$ が得られた時の予測分布 $P(\{ x_{N+t} \} | D)$ は，以下のようになる。

$\begin{align} P(\{ x_{N+t} \} | D) \\ &= \sum_{H} \iiint P(\{ x_{N+t} \} | \beta, H) P(w, \alpha, \beta, H | D) dw d\alpha d\beta \\ &= \sum_{H} \iiint P(\{ x_{N+t} \} | \beta, H) P(w |D, \alpha, \beta, H ) P(\alpha, \beta | D, H) P(H | D) dw d\alpha d\beta \\ \end{align}$

ややこしく見えるが，よくある予測分布の導出であり，パラメータ $H, w, \alpha, \beta$ を復活させて，離散パラメータは総和をとることで，連続パラメータは積分することで消去している。

グラフィカルモデルは，以下のようになる。

予測分布の近似

上記のように，予測分布はかなり複雑である。特に $P(w, \alpha, \beta, H | D)$ の計算が，一般に不可能に近いことが多い。そのため，以下のような近似を行なう。

$w$ : 事後分布最大値(MAP推定量)における2次近似を用いて積分する
$\alpha, \beta$ および $H$ : 事後分布最大値(MAP推定量)に固定する(そのため分布ではなくなる)。

パラメータ $w$ の２次近似について説明する。 $P(w |D, \alpha, \beta, H )$ を具体的に書き下すと，

$\begin{align} P(w |D, \alpha, \beta, H ) &= \frac{P(D | w, \beta, H ) P(w | \alpha, H) }{ P(D | \alpha, \beta, H)} \\ &= \frac{ \frac{1}{Z_W(\alpha) Z(\beta)} \exp (-M(w))} { \int \frac{1}{Z_W(\alpha) Z(\beta)} \exp (-M(w)) dw } \\ \end{align}$

ただし，

$\begin{align} M(w) &= \sum_{t=1}^{T-\tau} \left( -\frac{\beta}{2} (x_{t+\tau} - f(x_{t+\tau-1}, \cdots, x_t; w) ) \right) ^2 + \sum_{c=1}^{C} \frac{\alpha_c}{2} \lVert w_c \rVert ^2 \\ Z_W(\alpha) &= \prod_{c=1}^{C} Z_W(\alpha_c) , \quad Z(\beta) = Z_D(\beta)^{N-\tau} \\ \end{align}$

である。

$M(w)$ をMAP推定量 $w_{MP}$ の周りでTaylor展開すると，

$\begin{align} M(w) &\approx M(w_{MP}) + \frac{1}{2} (w - w_{MP})^T A (w - w_{MP}) \\ A &= \left. \frac{\partial ^2}{\partial w^2} M \right| _{w=w_{MP}} \\ \end{align}$

となる( $A$ は上記のようなヘッセ行列である)。

上式より，

$\begin{align} \exp(-M(w)) &\approx \exp(-M(w_{MP})) \exp \left( -\frac{1}{2} (w - w_{MP})^T A (w - w_{MP}) \right) \\ \end{align}$

さらに両辺を $w \in R^k$ で積分し，多変量正規分布の積分と正規化定数の関係を用いると，

$\begin{align} \int \exp(-M(w)) dw &\approx \exp(-M(w_{MP})) \int \exp \left( - \frac{1}{2} (w - w_{MP})^T A (w - w_{MP}) \right) dw \\ &= \exp(-M(w_{MP})) \times (2 \pi)^{k/2} \mathrm{det} (A^{-1})^{1/2} \\ \end{align}$

となる。これを用いると，

$\begin{align} P(w |D, \alpha, \beta, H ) &= \frac{ \frac{1}{Z_W(\alpha) Z(\beta)} \exp (-M(w))} { \int \frac{1}{Z_W(\alpha) Z(\beta)} \exp (-M(w)) dw } \\ \\ &\approx \frac{ \exp(-M(w_{MP})) \exp \left( -\frac{1}{2} (w - w_{MP})^T A (w - w_{MP}) \right) } { \exp(-M(w_{MP})) \times (2 \pi)^{k/2} \mathrm{det} (A^{-1})^{1/2} } \\ \\ &= \frac{1}{ (2 \pi)^{k/2} \mathrm{det} (A^{-1})^{1/2} } \exp \left( -\frac{1}{2} (w - w_{MP})^T A (w - w_{MP}) \right) \\ \end{align}$

となり， $P(w |D, \alpha, \beta, H )$ の２次近似が多変量正規分布として得られた*1。

最後に本節では， $\{ x_{N+t} \}$ の予測分布の近似についても言及している。予測問題においては，予測平均を用いることが多いため，予測分布を予測平均で近似する，としている。

第4章具体的問題

本節では，第3章までに定式化した予測手法を，以下の具体的問題に適用している。

ノイズを含むカオス的時系列予測(エノン系・レスラー系)
空調負荷予測

空調負荷予測については，空気調和・衛生工学会により実施された「熱負荷予測コンテスト」のデータを用いている。データセットは見つけられなかったが，コンテストの結果は以下に紹介されているものと考えられる。

http://www.tohoku-kuei.com/pdf/h28_gijutsu21.pdf

第5章ハミルトニアン・モンテカルロによるベイズ的学習と予測

第3章では，事後分布を近似式で表現していたが，本章ではMCMCによってサンプルを取得し，予測分布を求める方法について説明していた。

本書に関するまとめと感想

モデル選択も含めて階層ベイズモデルの枠組みで行なっていることが参考になった

本書は「階層ベイズモデル」がテーマであるが，非線形ダイナミカルシステム(非線形時系列モデル)をベイズ統計の枠組みで説明しきっている点が面白かった。

また事後分布についても，第3章では近似を用いて評価している点も参考になった。最近では，ベイジアンモデリングのソフトウェアが充実しているので，モデリングさえしっかりすれば実装はできるのであろうが，理論的な解析をしようとしたときにこのような近似は役立つのでは，と感じている。

実装にチャレンジしたい

今回のテーマは，ニューラルネットワークモデルに予測分布を与えるという内容である。これまで学んだ書籍は，基本的に線形モデルが多かったので，ベイズ統計と深層学習の融合分野やその実装については，須山敦志著「ベイズ深層学習」などを参考にしたい。

また実装についても，森賀新／木田悠歩／須山敦志著「Pythonではじめるベイズ機械学習入門」などを参考にしていきたい。

本記事を最後まで読んでくださり，どうもありがとうございました。

ランキング参加中

統計

*1:ただし，実際はヘッセ行列が出てくるので，パラメータ $w$ の次元数が大きくなりすぎると計算量は大きくなる

jiku log

データサイエンスの核心を掴む : 学びと発見の記録

「階層ベイズモデルとその周辺」を読む～第II部非線形ダイナミカルシステムの再構成と予測～