状態空間モデルの隠れ変数が連続的なモデルである線形動的システム(linear dynamical system, LDS)について紹介します。状態遷移確率、出力確率が条件付き線形ガウス分布で書けるという仮定が必要です。

具体的には以下の形でかけます。

$\begin{aligned} p(z_t | z_{t-1} ) &= \textit{N} ( z_{t} | A z_{t-1}, \Gamma)\\ p(x_{t} | z_{t}) &= \textit{N} ( x_{t} | C z_{t}, \Sigma)\\ p(z_1) &= \textit{N} (z_1 | \mu_0, P_0) \end{aligned}$

ここで、 $\textbf{X}, \textbf{Z}$ は観測変数系列、隠れ変数系列を表します。多くの場合では時系列データを考えるので、indexをtでおいて $\textbf{X} = ( x_1, \dots, x_t, \dots, x_T)$ とします。観測データの次元をN, 系列の長さをT, 隠れ変数の次元をMとします。 $x_t, z_t$ は連続値のベクトルになります。

グラフィカルモデルは以下のような形です。

f:id:endosan:20190730194224p:plain — 線形動的システムのグラフィカルモデル. z: 隠れ状態, x: 観測状態

線形動的システムにおいてもHMMと同じような問題が考えられます。

$p(X|\theta)$ を最大にするようなパラメーター $\theta = (A, \Gamma, C, \Sigma, \mu_0, P_0)$ を推定する
観測データを得た時に最適な隠れ状態の系列を知りたい

まずは２の問題について考えます。

最適な経路
- カルマンフィルター
  - 解釈
- カルマンスムーザ
  - 解釈
パラメーターの推定
- E step
- M step
  - パラメータの更新
補足
参考文献

最適な経路

2の問題については線形ガウス分布でモデル化されていることから、個々の時刻において最適な隠れ状態を取れば、最適な経路になることが言えます。そのためviterbiアルゴリズムのようなことはしなくてよくて、フォーワードバックワードアルゴリズムだけを計算すればよいです。

カルマンフィルター

カルマンフィルター(Kalman filter)として知られるフォーワードアルゴリズムについて取り上げます。ここで目標としているのは、観測データ $x_1, \dots, x_{t-1}$ を得た時に、次の時刻 $z_t, x_t$ の予測をすることです。

HMMでのスケーリングされたフォーワードアルゴリズムにしたがって進めます。観測データ $x_1, \dots, x_t$ を得た時に、隠れ状態 $z_t$ を得る確率について考えます。全ての条件付き確率がガウス分布なので次の分布はガウス分布でかけます。

$\begin{aligned} \hat{\alpha} (z_t) &:= p(z_t | x_1, \dots, x_t, \theta)\\ &= \textit{N} (z_t | \mu_t, V_t) \end{aligned}$

ここにおける $\mu_t, V_t$ はある定数です。

HMMでのフォーワードの式の積分版を考えると、

$\begin{aligned} c_t \hat{\alpha} (z_t) = p(x_t | z_t, \theta) \int \hat{\alpha} (z_{t-1}) p(z_t | z_{t-1}, \theta) d z_{t-1} \end{aligned}$

とかけます。ここで $c_t = p(x_t | x_{t-1}, \dots, x_1)$ のこと。これに具体的なガウス分布を代入してみると、( $\mu_t, V_t$ の具体的な形は今の段階では不明)

$\begin{aligned} c_t \textit{N} (z_t | \mu_t, V_t) &= \textit{N} ( x_{t} | C z_{t}, \Sigma) \int \textit{N} (z_{t-1} | \mu_{t-1}, V_{t-1}) \textit{N} ( z_{t} | A z_{t-1}, \Gamma) d z_{t-1}\\ c_t \textit{N} (z_t | \mu_t, V_t) &= \textit{N} ( x_{t} | C z_{t}, \Sigma) \textit{N} (z_t | A \mu_{t-1}, AV_{t-1}A^T + \Gamma)\\ c_t \textit{N} (z_t | \mu_t, V_t) &= \textit{N} ( x_{t} | C z_{t}, \Sigma) \textit{N} (z_t | A \mu_{t-1}, P_{t-1}) \end{aligned}$

where

$\begin{aligned} P_{t-1} = AV_{t-1}A^T + \Gamma \end{aligned}$

ガウス分布なので公式のようなものが利用できる

$\begin{aligned} c_t \textit{N} (z_t | \mu_t, V_t) = \textit{N} ( x_{t} | CA \mu_{t-1}, \Sigma+C P_{t-1} C^{T}) \textit{N} (z_{t} | A\mu_{t-1} + K_{t} (x_t - CA\mu_{t-1}), (I - K_t C) P_{t-1}) \end{aligned}$

ここで、Kalman gain matrix $K_t$ を次のように定義した。

$\begin{aligned} K_t = P_{t-1} C^{T} (C P_{t-1} C^{T} + \Sigma)^{-1} \end{aligned}$

よってパラメーターの時間発展として、次が得られます。

$\begin{aligned} \mu_t &= A\mu_{t-1} + K_{t} (x_t - CA\mu_{t-1})\\ V_t &= (I - K_t C) P_{t-1}\\ c_t &= \textit{N} ( x_{t} | CA \mu_{t-1}, \Sigma+C P_{t-1} C^{T}) \end{aligned}$

解釈

はじめにも言ったように、線形動的システムでは隠れ状態の最適な系列はその時刻ごとの最適な状態になるので、隠れ状態 $z_t$ の予測をするときには $A\mu_{t-1}$ を使います。そして観測状態 $x_t$ の予測は $CA\mu_{t-1}$ になります。これらは状態遷移確率、出力確率の仮定から従います。

次に、 $x_t$ を実際に観測すると、 $z_t$ は $x_t$ がgivenな状態に変わるので予測値 $\mu_t$ は

$\begin{aligned} \mu_t = A\mu_{t-1} + K_{t} (x_t - CA\mu_{t-1}) \end{aligned}$

のように修正が入って新しくなったと解釈できます。

カルマンスムーザ

カルマンスムーザー(Kalman smoother)はバックワードアルゴリズムに対応するものです。時刻Tまでの観測データを得ている時に隠れ状態 $z_t$ の推定をする時に使います。具体的には

$\begin{aligned} \gamma (z_t) = p(x_t | X, \theta) = \hat{\alpha}(z_t) \hat{\beta} (z_t) \end{aligned}$

を計算します。すでに $\hat{\alpha}(z_t)$ は求めているので、バックワードアルゴリズムから $\hat{\beta} (z_t)$ を求めます。また、フォーワードの計算の時と同じように、線形ガウスで全ての分布が書かれていることから

$\begin{aligned} \gamma (z_t) = \textit{N} (z_t | \hat{\mu}_t, \hat{V}_t ) \end{aligned}$

と置けることがわかります。

バックワードのスケーリングされた時間発展の式は以下のようになります。（HMMで求めたものを隠れ状態については積分形にしただけです）

$\begin{aligned} c_{t+1} \hat{\beta} (z_t) = \int \hat{\beta} (z_{t+1}) p(x_{t+1} | z_{t+1}) p(z_{t+1} | z_t) d z_{t+1} \end{aligned}$

両辺に $\hat{\alpha}(z_t)$ をかけると、

$\begin{aligned} c_{t+1} \hat{\alpha}(z_t)\hat{\beta} (z_t) &= \hat{\alpha}(z_t) \int \hat{\beta} (z_{t+1}) p(x_{t+1} | z_{t+1}) p(z_{t+1} | z_t) d z_{t+1}\\ c_{t+1} \gamma (z_t) &= \hat{\alpha}(z_t) \int \hat{\beta} (z_{t+1}) \textit{N} ( x_{t+1} | C z_{t+1}, \Sigma) \textit{N} ( z_{t+1} | A z_{t}, \Gamma) d z_{t+1}\\ \end{aligned}$

ここの過程はあとで書きます。

最終的にパラメーターの時間発展として、次が得られます。

$\begin{aligned} \hat{\mu}_t &= \mu_{t} + J_{t} (\hat{\mu}_{t+1} - A\mu_{t})\\ \hat{V}_t &= V_t + J_t (\hat{V}_{t+1} - P_t) J_t^{T} \end{aligned}$

解釈

全時刻の観測データを使うので、より正確な隠れ状態の推定ができていると考えることができます。隠れ状態の平均 $\hat{\mu}_t$ は、時刻1からtまで得た時の推定値 $\mu_t$ に修正 $J_{t} (\hat{\mu}_{t+1} - A\mu_{t})$ を加えたと考えることができます。

パラメーターの推定

観測データを時刻Tまでえたらパラメータの推定ができます。厳密な解析解があるわけではないので、EMアルゴリズムによってパラメーターを求めます。

完全データの対数尤度関数の値が収束したら終了です。現実の利用では、パラメーターの変化がある程度小さくなったら打ち切りにすることが多いかと思います。

E step

KL divergenceを０にする隠れ変数 $Z$ の事後分布は現在の尤度関数 $p$ の形からわかります。

$\begin{aligned} q(Z) = p(Z | X, \theta_{old}) \end{aligned}$

M stepで必要な期待値の計算をします。

$\begin{aligned} E[ z_t] &= \hat{\mu}_t\\ E[z_t z_{t-1}^{T}] &= \hat{V}_t J_{t-1}^{T} + \hat{\mu}_t \hat{\mu}_{t-1} \end{aligned}$

M step

Q関数をパラメーター $\theta$ で微分することにより求められる。まずはQ関数の具体的形をみる。そのために完全データの対数尤度をかくと以下の形になる。

$\begin{aligned} \log p(X, Z | \theta) &= \log p(z_1 | \mu_0, P_0) + \sum_{t=2}^{T} \log p (z_t | z_{t-1}, A, \Gamma) + \sum_{t=1}^{T} \log p(x_t | z_t, C, \Sigma) \end{aligned}$

するとQ関数は

$\begin{aligned} Q ( \theta, \theta_{old}) & = E_{Z|\theta_{old}} \left[ \log p(X, Z | \theta)\right] \end{aligned}$

となる。具体的な分布を代入すると以下の形になる。

$\begin{aligned} Q ( \theta, \theta_{old}) =& -\frac{1}{2} \log |P_0| - E_{Z | \theta_{old}} [ \frac{1}{2} ( z_1 - \mu_0)^T P_{0}^{-1}(z_1 - \mu_0) ]\\ & -\frac{T-1}{2} \log |\Gamma| - E_{Z | \theta_{old}}[ \frac{1}{2} \sum_{t=2}^{T} z_{t} - A z_{t-1})^T \Gamma^{-1} (z_{t} - A z_{t-1}) ] \\ & - \frac{T}{2} \log |\Sigma| - E_{Z | \theta_{old}} [ \frac{1}{2} \sum_{t=1}^{T} (x_{t} - Cz_{t})^T \Sigma^{-1} (x_{t} - Cz_{t}) ] + \textrm{const.} \end{aligned}$

パラメータの更新

Q関数をパラメーター $\theta = (A, \Gamma, C, \Sigma, \mu_0, P_0)$ について微分をする。

$\begin{aligned} \frac{\partial Q ( \theta, \theta_{old})}{\partial \mu_0} = 0 &\rightarrow \mu_0 = E_{Z | \theta_{old}} [ z_1 ] \\ \frac{\partial Q ( \theta, \theta_{old})}{\partial P_0} = 0 &\rightarrow P_0 = E_{Z | \theta_{old}} [ z_1 z_1^T] - E_{Z | \theta_{old}} [ z_1 ] E_{Z | \theta_{old}} [ z_1^T ]\\ \frac{\partial Q ( \theta, \theta_{old})}{\partial A} = 0 &\rightarrow A = \left( \sum_{t=2}^{T} E_{Z | \theta_{old}} [ z_t z_{t-1}^T] \right) \left( \sum_{t=2}^{T} E_{Z | \theta_{old}} [ z_{t-1} z_{t-1}^T] \right)^{-1}\\ \frac{\partial Q ( \theta, \theta_{old})}{\partial \Gamma} = 0 &\rightarrow \Gamma = \frac{1}{T-1} \sum_{t=2}^{T} \left\{ E_{Z | \theta_{old}} [z_{t} z_{t}^{T}] - A E_{Z | \theta_{old}} [z_{t-1} z_{t}^{T}] - E_{Z | \theta_{old}} [z_{t} z_{t-1}^{T}] A^{T} + A E_{Z | \theta_{old}} [z_{t-1} z_{t-1}^{T}] A^{T} \right\}\\ \frac{\partial Q ( \theta, \theta_{old})}{\partial C} = 0 &\rightarrow C = \left( \sum_{t=1}^{T} x_{t} E_{Z | \theta_{old}} [z_t^T]\right) \left( \sum_{t=1}^{T} E_{Z | \theta_{old}} [z_{t} z_{t}^{T}] \right)^{-1}\\ \frac{\partial Q ( \theta, \theta_{old})}{\partial \Sigma} = 0 &\rightarrow \Sigma = \frac{1}{T} \sum_{t=1}^{T} \left\{ x_{t} x_{t}^{T} - C E_{Z | \theta_{old}} [z_{t} ] x_{t}^{T} - x_{t} E_{Z | \theta_{old}}[z_{t}^{T}] C^{T} + C E_{Z | \theta_{old}}[z_{t} z_{t}^{T}] C^{T} \right\} \end{aligned}$

補足

ガウス分布の便利な計算について

$\begin{aligned} p(x) &= \textit{N} (x | \mu, \Lambda^{-1})\\ p(y|x) &= \textit{N} (y | Ax + b, L^{-1}) \end{aligned}$

ならば、

$\begin{aligned} p(y) &= \textit{N} (y | A\mu+b, \Lambda^{-1} + A \Lambda^{-1} A^{T})\\ p(y|x) &= \textit{N} (x | \Sigma \{ A^{T} L ( y-b)+ A\mu \}, \Sigma) \end{aligned}$

where

$\begin{aligned} \Sigma = ( \Lambda + A^{T} L A)^{-1} \end{aligned}$

行列計算について

$\begin{aligned} (P^{-1} + B^{T} R^{-1} B)^{-1} B^{T} R^{-1} &= P B^{T} (BPB^{T} + R)^{-1}\\ (A+BD^{-1}C)^{-1} &= A^{-1} - A^{-1} B (D+CA^{-1}B)^{-1} CA^{-1} \end{aligned}$