カルマンフィルタ リング問題に代表されるような、線形動的システム(Linear Dynamical System)のパラメータをEMアルゴリズムで推定してみます。
やってる内容そのものは、「Parameter Estimation for Linear Dynamical Systems」の実装になります。この論文は、25年前に現Google BrainのテックリードであるGhahramani氏とDeepLearningの親と呼ばれるHintonがが1996年に出した論文です。
内容自体は基本的で、PRMLに掲載されているような機械学習の基本的な内容だと思いますが、抽象的なとこもあり理解も大変なので、この論文の内容を解説したのちに、Pythonで実装してみたいと思います。
- 論文「Parameter Estimation for Linear Dynamical Systems」の解説
- 線形動的モデルのパラメータ推定をEMアルゴリズムで行う
- PythonでEMアルゴリズムの実装と線形動的モデルの実装を行う
今回考える線形動的モデル
論文の内容に即して解説していきます。
今回考えるモデルは、下記のような一般的な線形動的システムです。
状態方程式と観測方程式
\begin{equation}
\begin{split}
\bm{x}_{t} &= A\bm{x}_{t-1} + \bm{w}_t \\
\bm{y}_t &= C\bm{x}_t + \bm{v_t}
\end{split}
\end{equation}ここで、$w_t$、$v_t$は平均0、分散共分散がそれぞれ$R , Q$のガウス分布に従うホワイトノイズであるとします。
\begin{equation}
\begin{split}
\bm{w}_t &\sim \mathcal{N}(0, R) \\
\bm{v}_t &\sim \mathcal{N}(0, Q)
\end{split}
\end{equation}時刻$t$によって分散共分散が変わるようなモデルでも対応できますが、今回は分散共分散は時刻によらず一定のモデルとしています。
また、状態$\bm{x}$の初期値$x_1$は、下記のように設定しています。
\begin{equation}
\begin{split}
\bm{x}_1 &\sim \mathcal{N}(\bm{\pi}_1, V_1) \\
\end{split}
\end{equation}(3)式によって、状態の初期値がガウス分布に従うことをモデルに導入しています。(1)式の状態方程式と観測方程式は、すべて線形変換によって構成されているので、ガウス分布の線型性の性質から、観測値も状態列もすべてガウス分布に従うことになります。
今回扱う、一般的な線形動的モデルに登場する全てのモデルパラメータ\bm{\theta}は、下記のようになります。
\begin{equation}
\bm{\theta} = \{ A, C, R, Q, \bm{π_1}, V_1\}
\end{equation}これらのパラメータを全てEMアルゴリズムで学習していくことが、本記事の目的となります。
線形動的モデルの確率的な表現
ここまでで、線形動的モデルの一般的な表現について見てきまた。続いて、(1)~(3)のモデルを確率的に表現していきましょう。
まず、(1)式の2番目に登場する観測方程式を確率的な表現にしてみましょう。
観測方程式の表現から、条件付き確率分布 $p(y_t| x_t)$ を次のように表現することができます。
\begin{equation}
\begin{split}
p(y_t | x_t) &= exp
\biggr\{
- \frac{1}{2}(y_t - Cx_t)^TR^{-1}(y_t-Cx_t)
\biggr\}
{(2\pi)}^{-\frac{p}{2}} {|R|}^{- \frac{1}{2}} \\
&= - \frac{1} {(\sqrt{2\pi})^p \sqrt{|R|}}
exp
\biggr\{
- \frac{1}{2}(y_t - Cx_t)^TR^{-1}(y_t-Cx_t)
\biggr\}
\end{split}
\end{equation}この確率は、時刻$t$の状態$\bm{x}_t$がわかっている時、$\bm{y}_t$に遷移する確率を示しています。
この式変形は、ガウス分布の線形変換は、ガウス分布になることを利用しています。
続いて、観測行列の表現から、条件付き確率分布$p(x_t| x_{t-1})$を次のように確率的に表現することができます。
\begin{equation}
\begin{split}
p(x_t | x_{t-1}) &= exp
\biggr\{
- \frac{1}{2}(x_t - Ax_{t-1})^TQ^{-1}(x_t- Ax_{t-1})
\biggr\}
{(2\pi)}^{-\frac{k}{2}} {|R|}^{- \frac{1}{2}} \\
&= - \frac{1} {(\sqrt{2\pi})^k \sqrt{|Q|}}
exp
\biggr\{
- \frac{1}{2}(x_t - Ax_{t-1})^TQ^{-1}(y_t-Ax_{t-1})
\biggr\}
\end{split}
\end{equation}ここでも、ガウス分布の線形変換が、ガウス分布になることを利用しました。
続いて、(3)式で示されている初期値についても確率的な表現に直します。
線形動的モデルではガウシアンを仮定しているので、初期値$x_1$の平均を$\pi_1$分散共分散行列$V_1$を用いて、初期状態は次のような確率で表現することができます。
状態の初期値$x_1$の平均を$\pi_1$、分散共分散行列である$V_1$とすると、初期状態$x_1$の確率密度関数は次のように表現することができる。
\begin{equation}
\begin{split}
p(x_1) &= exp
\biggr\{
- \frac{1}{2}(x_1 - \pi_1)^TV^{-1}(x_1- \pi_1)
\biggr\}
{(2\pi)}^{-\frac{k}{2}} {|V_1|}^{- \frac{1}{2}} \\
&= - \frac{1} {(\sqrt{2\pi})^k \sqrt{|V_1|}}
exp
\biggr\{
- \frac{1}{2}(x_1 - \pi_{1})^TV_1^{-1}(x_1-\pi_{1})
\biggr\}
\end{split}
\end{equation}ここまでで、やっと確率的な表現をする準備まですることができました。
続いて、EMアルゴリズムの内容に入っていきます。
線形動的モデルにおけるEMアルゴリズム
それでは、まず線形動的システムにおけるEMアルゴリズムの内容に入っていきましょう。
線形動的システムにおいては、観測値$\bm{Y}$が得られている中で、(4)式で示されているような全パラーメータ$\bm{\theta}$を学習していくことになります。
EMアルゴリズムの導入の内容にもなりますが、得られたデータ$\bm{Y}$から、尤度関数$p(\bm{Y} | \bm{\theta})$を用いてパラメータ$\bm{\theta}$の最適化を行うことは、解析的に行うことは難しいので、EMアルゴリズムでは、潜在変数と観測値の同時分布で尤度関数を表現し、その同時分布を扱うことで、パラメータの最適化を行います。
尤度関数$p(\bm{Y} | \bm{\theta})$を潜在変数と観測値の同時分布で表現するとこのようになります。
\begin{equation}
\begin{split}
p(\bm{Y} | \bm{\theta}) = \int_x p(\bm{X}, \bm{Y} | \bm{\theta}) dx
\end{split}
\end{equation}(8)式の右辺を最適化することで、最適なパラメータを計算していきます。
EMアルゴリズムにおける具体的な計算
EMアルゴリズムでは、潜在変数における完全データの対数尤度関数の期待値を計算するEステップと、Eステップによって得られた$\mathcal{Q} (\bm{\theta}, \bm{\theta^{old}})$関数を、各パラメータ$\bm{\theta}$で最適化するMステップに分かれます。
EMアルゴリズムにおけるEステップ
Eステップでは、完全データ対数尤度の期待値を計算することがゴールです。具体的には(9)式を導出することが重要です。
完全データ対数尤度を、$ln p(\bm{X}, \bm{Y} | \bm{\theta})$とした時、EMアルゴリズムによって最適化する、完全データ対数尤度の期待値$\mathcal{L}( \bm{\theta}, \bm{\theta}^{old})$は下記のように表現されます。
\begin{equation}
\begin{split}
\mathcal{Q} (\bm{\theta}, \bm{\theta^{old}}) =
\int_X p(\bm{X} | \bm{Y}, \bm{\theta^{old}})ln p(\bm{Y}, \bm{X} | \bm{\theta}) dx
\end{split}
\end{equation}Eステップでやっている 内容は、まさにカルマンフィルタリングとカルマンスムージング(Rauch Tung Striebel RTSスムージング)です。これらの内容を見ていきましょう。
Eステップの計算の準備
まず、Eステップの準備をしていきます。今、簡便さのためパラメータ$\bm{\theta^{old}}$と観測値$\bm{Y}$の条件のもとでの、時刻$t$の状態の確率を次のように定義します。
状態の事後分布はガウス分布で表現できるので、次の(10)式のように定義できる。
\begin{equation}
\begin{split}
\alpha(\bm{x_n}) &= p(\bm{x_n} | \bm{y_1}, \bm{y_2}, \dots, \bm{y_n}, \bm{\theta}^{old}) \\
&= \mathcal{N} (\bm{x_n} | \mu_n, V_n)
\end{split}
\end{equation}時刻$t$の状態の事後分布を、$\alpha(\bm{x_n})$のように定義しました。
$\alpha(\bm{x_n})$は、ガウス分布に従うため、そのパラメータ$\mu_n, V_n$とおると、これらのパラメータ$\mu_n, V_n$は、次のようなアルゴリズムで再帰的に導出できることが知られています。
\begin{equation}
\begin{split}
\mu_n = A^{old} \mu_{n-1} + K_{n-1} (\bm{x_n} - C^{old}A^{old}\mu_{n-1})
\end{split}
\end{equation}\begin{equation}
\begin{split}
V_n = (I - K_{n-1} C^{old})P_{n-1}
\end{split}
\end{equation}ここで、$P_n$、$K_n$の定義はこのようになっている。
\begin{equation}
\begin{split}
P_n = A^{old}V_n(A^{old})^T + Q
\end{split}
\end{equation}\begin{equation}
\begin{split}
K_n = P_n(C^{old})^T (C^{old}P_n (C^{old})^T + R^{old})^{-1}
\end{split}
\end{equation}(11)、(12)で再帰的に計算がされているが、初期値についてはこのようになっている。
\begin{equation}
\begin{split}
\mu_0 = \pi_1^{old} + K_0(x_0 - C^{old} \pi_1^{old})
\end{split}
\end{equation}\begin{equation}
\begin{split}
V_0 = (I - )
\end{split}
\end{equation}また、線形ガウシアンモデルではマルコフ性が仮定されているので、全状態変数$\bm{x}$と全観測値$\bm{y}$の同時確率分布は次のように定義されます。
\begin{equation}
p(\bm{x}, \bm{y}) = p(x_1) \prod_{t=2}^{T}p(x_t| x_{t-1}) \prod_{t=1}^T(y_t | x_t)
\end{equation}これは、条件付き確率分布の連鎖律(Chain Rule)の公式を使うことで、同時分布を条件付き確率を用いてこのように書き下すことがでいます。条件付き確率分布の連鎖律に関してはこちらの記事で解説しています。
同時確率の対数を考える
(5)式で、全部の状態、全部の観測値の同時確率分布を得ることができました。続いて、この同時確率分布は積の形をして分かりにくいので、この対数をとります。
元の論文の(7)式で表現されているのですが、ここはかなり式変形がカオスになるので、ゆっくりと順番を追ってやっていきます。
(5)式で表現した同時確率分布の対数$ln p(\bm{x}, \bm{y})$は次のようになる。
\begin{equation}
\begin{split}
ln p(\bm{x}, \bm{y}) &= lnp(x_1) + ln
\biggr\{
\prod_{t=2}^{T}p(x_t| x_{t-1})
\biggr\} +
ln \biggr\{
\prod_{t=1}^T p(y_t | x_t)
\biggr\} \\
&=
lnp(x_1) +
\sum_{t=2}^{T}
\biggr(
lnp(x_t| x_{t-1})
\biggr) +
\sum_{t=1}^{T}
\biggr(
ln p(y_t | x_t)
\biggr) \\
&=
- \sum_{t=1}^{T}
\biggr(
- \frac{1}{2}(y_t - Cx_t)^TR^{-1}(y_t-Cx_t)
\biggr)
-\frac{T}{2}log|R| \\
&\space\space
- \sum_{t=2}^{T}
\biggr(
- \frac{1}{2}(x_t - Ax_{t-1})^TQ^{-1}(x_t-Ax_{t-1})
\biggr)
-\frac{T-1}{2}log|Q| \\
&\space\space
- \frac{1}{2}(x_1 - \pi_{1})^TV_1^{-1}(x_1-\pi_{1})
- \frac{1}{2} log|V_1| - \frac{T(p+k)}{2}log2\pi
\end{split}
\end{equation}これが、論文中の(7)式と対応します。EMアルゴリズムの教科書に登場するような表記で記載するとこのようになります。
EMアルゴリズムにおけるMステップ
Mステップにおいては、Eステップで求めた事後分布を利用できるので、完全データ対数尤度の期待値を最大化しするような、パラメータ$\bm{\theta}$を求めます。
線形動的モデルなので、パラメータ$A$, $C$は線形変換を行う係数ベクトルとなり、$w_t, v_t$は、それぞれ、平均が0、分散が$R, Q$のガウス分布に従います。
ここでモデルパラメータを$\bm{\theta}$とすると、今回のパラメータは下記のようなものがあります。
EMアルゴリズムにおける、Mステップでは、尤度関数$\mathcal{L}$に対し、(3)式で示したモデルパラメータで偏微分し、1ステップ毎の最尤推定値を導出していきます。
この際、$x_{t}$が潜在変数となるため、このモデルにおいては、$x_{t}$を固定し、この値を固定した上で、(3)式の対数尤度関数をパラメータで偏微分をし、パラメータの更新を行います。