隠れマルコフモデル
この記事は、須山敦志著の「ベイズ推論による機械学習入門」を参考に作成した。
いつ見返しても理解が滞らないよう、数式の扱いがくどい仕様になっている。
隠れマルコフモデルは、隠れ変数(潜在変数)にマルコフ性を仮定したモデルであり、時系列データの分析に用いられる結構有名なモデルらしい。
ここでは、一つ前の時刻のデータのみに依存する1次マルコフ連鎖を扱う。
混合モデルでは潜在変数 $ S $ の発生には条件付き独立を仮定していた。
$$
p(\textbf{S}|\textbf{π}) = \prod_{n=1}^{N} p(\textbf{s}_n|\textbf{π})
$$
隠れマルコフモデルでは隣り合う潜在変数 $ \textbf{s}_n, \textbf{s}_{n-1}$ に依存性を設ける。
$$
p(\textbf{S}|\textbf{π}, \textbf{A})
= p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})
$$
隠れマルコフモデルでは、潜在変数 $ \textbf{S} $ のことを 状態系列 と呼ぶ。
ここで、$\textbf{A} \in \mathbb{R}^{K × K} $ は遷移確率行列と呼ばれ、$A_{i,j}$は状態$j$から状態$i$に遷移する確率を表す。ここで、$K$は状態数(クラスタ数) $K$ である。
たとえば、$K = 3$ の $A$ の一例は以下である。
\textbf{A}=
\begin{bmatrix}
0.5 & 0.0 & 0.1 \\
0.2 & 0.8 & 0.4 \\
0.3 & 0.2 & 0.5 \\
\end{bmatrix}
遷移確率行列は、列ごとにみると分かりやすい。
$$
\sum_{j=1}^{K} \textbf{A}_{j, i}
$$
であることが、遷移確率行列の意味からも、上記の$\text{A}$の例からもわかる。
$\textbf{s}_1$ は $\textbf{π}$ をパラメータとするカテゴリ分布から発生すると考える。
$\textbf{π}$ は $\textbf{s}_1$ のみに影響を与えるため、初期確率と呼ばれる。
$$
p(\textbf{s}_1|\textbf{π}) = Cat(\textbf{s}_1|\textbf{π})
$$
また、カテゴリ分布のパラメータ $\textbf{π}$ の事前分布には、カテゴリ分布の共役事前分布であるディリクレ分布を採用する。
$$
p(\textbf{π}) = Dir(\textbf{π}|\textbf{α})
$$
ここで、隠れマルコフモデルの潜在変数と遷移確率行列について改めて考える。
$$
p(\textbf{S}|\textbf{π}, \textbf{A})
= p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})
$$
の $p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})$ について。
\begin{eqnarray*}
p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})
&=& \prod_{i=1}^{K} Cat(\textbf{s}\_n | \textbf{A}_{:,i})^{\textbf{s}_{n-1}, i} \\
&=& \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j,i}^{\textbf{s}_{n,j} ~ \textbf{s}_{n-1,i}} \\
\end{eqnarray*}
となる。$\textbf{A}_{:,i}$ は、$\textbf{A}$ の $i$ 行目のベクトルを表す。
ここは用心して、先ほどの例
\textbf{A}=
\begin{bmatrix}
0.5 & 0.0 & 0.1 \\
0.2 & 0.8 & 0.4 \\
0.3 & 0.2 & 0.5 \\
\end{bmatrix}
で具体的に計算してみる。
例えば、$\textbf{s}_{1} = (0,1,0)^T$ のときに $\textbf{s}_2 = (0,0,1)^T $ となる確率を求める。
\begin{eqnarray*}
p(\textbf{s}_2 | \textbf{s}_1, \textbf{A})
&=& \prod_{i=1}^{3} Cat((0,0,1)^T | \textbf{A}_{:,i})^{\textbf{s}_{n-1, i}} \\
&=& Cat((0,0,1)^T | (A_{1,1}, A_{2,1}, A_{3,1}))^0・Cat((0,0,1)^T | (A_{1,2}, A_{2,2}, A_{3,2}))^1・Cat((0,0,1)^T | (A_{1,3}, A_{2,3}, A_{3,3}))^0 \\
&=& Cat((0,0,1)^T | (A_{1,2}, A_{2,2}, A_{3,2}))^1 \\
&=& A_{3,2}
\end{eqnarray*}
なるほど、$\textbf{s}_n$ の要素が $0$ か $1$ であれば、
p(\textbf{s}_{n,j} = 1 | \textbf{s}_{n-1,i}, \textbf{A}) = A_{i,j}
となるわけだ。
事後予測として $\textbf{s}_{n,k}$ の要素が $0$ から $1$ の間の実数値になると、$ p(\textbf{s}_2 | \textbf{s}_1, \textbf{A}) $ はもう少し複雑になるのかな。
$\textbf{A}$ の $i$ 番目の列 $\textbf{A}_{:,i}$ は状態 $i$ から状態 $j (j= 1, 2, .. , K)$ を選択するカテゴリ分布であり、$\textbf{A}_{:,i}$ の事前分布は、カテゴリ分布の共役事前分布であるディリクレ分布を採用する。
p(\textbf{A}_{:,i}) = Dir(\textbf{A}_{:,i} | \textbf{β}_{:,i})
$\textbf{β}$ は $\textbf{A}$ と同じ $K×K$ の行列である。
隠れマルコフモデルの同時分布を構成する準備はできた。
モデルのパラメータを $\textbf{Θ}$ とすれば、同時分布は以下となる。
\begin{eqnarray*}
p(\textbf{X}, \textbf{S}, \textbf{Θ}, \textbf{π}, \textbf{A})
&=& p(\textbf{X} | \textbf{S}, \textbf{Θ}) p(\textbf{S} | \textbf{π}, \textbf{A}) p(\textbf{Θ}) p(\textbf{π}) p(\textbf{A}) \\
&=& p(\textbf{Θ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{x}_1|\textbf{s}_1, \textbf{Φ}) p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{x}_n | \textbf{s}_n, \textbf{Θ}) p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) \\
\end{eqnarray*}
今回は、ポアソン分布を仮定してモデルをつくる。
p(x_n | \textbf{s}_n, \textbf{λ}) = \prod_{k=1}^{K} Poi(x_n | λ_k)^{s_{n,k}}
パラメータ $ λ_k $ の事前分布には、ポアソン分布の共役事前分布であるガンマ分布を採用する。
p(λ_k) = Gam(λ_k | a, b)
今回使う分布の紹介
カテゴリ分布
ベルヌーイ分布の$K$次元への拡張。$π_k \in ( 0,1 )$、$\sum_{k=1}^{K} π_k = 1$、$s_k \in { 0,1 }$、$\sum_{k=1}^{K} s_k = 1$ である。
$$
Cat(\textbf{s} | \textbf{π})
= \prod_{k=1}^{K} π_k^{s_k}
$$
$$
\ln Cat(\textbf{s} | \textbf{π})
= \sum_{k=1}^{K} s_k \ln π_k
$$
ポアソン分布
非負の整数 $x$ を生成するための確率分布。$λ \in \mathbb{R}^+$ である。
$$
Poi(x|λ)
= \frac{λ^x}{x!} e^{-λ}
$$
$$
\ln Poi(x|λ)
= x \ln λ - \ln x! - λ
$$
ディリクレ分布
ハイパーパラメータ $ \textbf{α} \in \mathbb{R}^{+K} $ に基づき、$K$ 次元ベクトル $ \textbf{π} $ を生成する。
$$
Dir(\textbf{π}|\textbf{α})
= C_D(\textbf{α}) \prod_{k=1}^{K} π_k^{α_k - 1}
$$
$$
\ln Dir(\textbf{π}|\textbf{α})
= \sum_{k=1}^{K} (α_k - 1) \ln π_k + \ln C_D(\textbf{α})
$$
ただし、正規化項は以下である。
$$
C_D(\textbf{α}) = \frac{Γ(\sum_{k=1}^{K}α_k)}{\prod_{k=1}^{K} Γ(α_k)}
$$
ガンマ分布
正の実数 $λ \in \mathbb{R}^+$ を生成する。ハイパーパラメータは $a,b \in \mathbb{R}^+$ である。
$$
Gam(λ|a,b) = C_G(a,b) λ^{a-1} e^{-b λ}
$$
$$
\ln Gam(λ|a,b) = (a-1) \ln λ - b λ + \ln C_G(a,b)
$$
ただし、正規化項は以下である。
$$
C_G(a,b) = \frac{b^a}{Γ(a)}
$$
変分推論
変分推論の手法を導入する。
確率分布 $p(z_1,z_2,z_3)$ を $q(z_1)q(z_2)q(z_3)$ で近似することを考える。
分解することで、$z_i$ と $z_j$ $(i≠j)$ の独立性が仮定され、相関をとらえられなくなる代わりに、より単純な式で表せたり、計算量の削減につながる。
$q(z_1)$ の近似の仕方は以下である。
$$
q_{opt.}(z_1) = \underset{q(z_1)}{argmin} KL[q(z_1)q(z_2)q(z_3) || p(z_1, z_2, z_3)]
$$
分解近似した分布と元の分布のKLダイバージェンスを $q(z_1)$ について最小化する。
以下、分布 q(z_m) による期待値の計算を $<・>_{q(z_m)}$ と表記する。
\begin{eqnarray*}
KL[q(z_1)q(z_2)q(z_3) || p(z_1, z_2, z_3)]
&=& \int q(z_1)q(z_2)q(z_3) \ln \frac{q(z_1)q(z_2)q(z_3)}{p(z_1, z_2, z_3)} dz_1 \\
&=& <\ln \frac{q(z_1)q(z_2)q(z_3)}{p(z_1, z_2, z_3)}>_{q(z_1)q(z_2)q(z_3)} \\
&=& <\ln q(z_1)>_{q(z_1)} + <\ln q(z_2)>_{q(z_2)} + <\ln q(z_3)>_{q(z_3)} - <\ln p(z_1, z_2, z_3)>_{q(z_1)q(z_2)q(z_3)} \\
&=& <\ln q(z_1)>_{q(z_1)} - <<\ln p(z_1, z_2, z_3)>_{q(z_2)q(z_3)}>_{q(z_1)} + const \\
&=& <\ln q(z_1)>_{q(z_1)} - <\ln exp\{<\ln p(z_1, z_2, z_3)>_{q(z_2)q(z_3)}\}>_{q(z_1)} + const \\
&=& <\ln \frac{q(z_1)}{exp\{<\ln p(z_1, z_2, z_3)>_{q(z_2)q(z_3)}\}}>_{q(z_1)} + const \\
&=& KL[q(z_1) || exp\{<\ln p(z_1, z_2, z_3)>_{q(z_2)q(z_3)}\}] + const \\
\end{eqnarray*}
$q(z_2)$ と $q(z_3)$ が固定されたとすると、$KL[q(z_1)q(z_2)q(z_3) || p(z_1, z_2, z_3)]$ を最小にする $q(z_1)$ は次の式から得られる。
$$
\ln q(z_1) = <\ln p(z_1, z_2, z_3)>_{q(z_2)q(z_3)} + const
$$
実際に使う時を考え、パラメータ $z_m (m=1,2,..,M)$ と観測データ $D$ についての近似について考える。
$M$ 個あるパラメータ $z_m$ のうち、$z_i$ のみを除いた集合を、$Z_{＼i}$ と表記する。
このモデルの事後分布 $p(z_1, z_2,.., z_M|D)$ を $q(z_1)q(z_2)･･･q(z_M)$ で近似する。
\begin{eqnarray*}
\ln q(z_i)
&=& <\ln p(z1, .., z_M | D)>_{q(Z_{＼i})} + const \\
&=& <\ln \frac{p(z1, .., z_M, D)}{p(D)}>_{q(Z_{＼i})} + const \\
&=& <\ln p(z1, .., z_M, D)>_{q(Z_{＼i})} + const \\
\end{eqnarray*}
上記のようになり、同時分布 $ p(z1, .., z_M, D) $ を用いて $q(z_i)$ を計算できる。
完全分解変分推論
パラメータの事後分布を分解近似する。
$$
p(\textbf{S}, \textbf{λ}, \textbf{A}, \textbf{π} | X) ≈ { \prod_{n=1}^{N} q(\textbf{s}_n) } q(\textbf{λ}, \textbf{A}, \textbf{π})
$$
上式では時系列方向の分解を仮定しており、このような近似を完全分解変分推論と呼ぶ。
同時分布を今一度ここで。
\begin{eqnarray*}
p(\textbf{X}, \textbf{S}, \textbf{λ}, \textbf{π}, \textbf{A})
&=& p(\textbf{X} | \textbf{S}, \textbf{λ}) p(\textbf{S} | \textbf{π}, \textbf{A}) p(\textbf{λ}) p(\textbf{π}) p(\textbf{A}) \\
&=& p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{S} | \textbf{π}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ}) \\
&=& p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ}) \\
\end{eqnarray*}
さっそく、変分推論してみよう。
\begin{eqnarray*}
\ln q(\textbf{λ}, \textbf{π}, \textbf{A})
&=& <\ln p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{S} | \textbf{π}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{S})} + const \\
&=& <\ln p(\textbf{λ}) + \ln p(\textbf{π}) + \ln p(\textbf{A}) + \ln p(\textbf{S} | \textbf{π}, \textbf{A}) + \sum_{n=1}^{N} \ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{S})} + const \\
&=& \ln p(\textbf{λ}) + \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{S} | \textbf{π}, \textbf{A})>_{q(S)} + \sum_{n=1}^{N} <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{s})_n} + const \\
\end{eqnarray*}
とりあえずここまで。
次は $ q(\textbf{s}_i) $ を変分推論する。ここで、少し注意。たいていの $ \textbf{s}_i $ は $ \textbf{s}_{i-1} $ に依存性をもち、$ \textbf{s}_{i+1} $ に依存性をもたせる。
たいていの場合でないのは、インデックスが端っこのものだ。$ \textbf{s}_1 $ は事前確率 $\textbf{π}$ に依存性をもち、$ \textbf{s}_2 $ に依存性をもたせる。$ \textbf{s}_N $ は $ \textbf{s}_{N-1} $ に依存性をもつのみ。
ということで、$i=1, i=n(n=2,..,N-1), i=N$三つのパターンに分けて変分推論を行う。
以下、着目するパラメータに関係ない項はどんどん $const$ に吸収させちゃおう。
$i=1$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_1)
&=& <\ln p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼1})} + const \\
&=& <\ln p(\textbf{s}_1 | \textbf{π}) + \ln p(\textbf{s}_2 | \textbf{s}_1, \textbf{A}) + \ln p(x_1 | \textbf{s}_1, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼1})} + const \\
&=& <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{π})} + <\ln p(\textbf{s}_2 | \textbf{s}_1, \textbf{A})>_{q(\textbf{s}_2)q(\textbf{A})} + <\ln p(x_1 | \textbf{s}_1, \textbf{λ})>_{q(\textbf{λ})} + const \\
\end{eqnarray*}
$i=n~(n=2,..,N-1)$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_n)
&=& <\ln p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼n})} + const \\
&=& <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) + \ln p(\textbf{s}_{n+1} | \textbf{s}_n, \textbf{A}) + \ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼n})} + const \\
&=& <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n-1})} + <\ln p(\textbf{s}_{n+1} | \textbf{s}_n, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n+1})} + <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} + const \\
\end{eqnarray*}
$i=N$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_N)
&=& <\ln p(\textbf{λ}) p(\textbf{π}) p(\textbf{A})・p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) \prod_{n=1}^{N} p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼n})} + const \\
&=& <\ln p(\textbf{s}_N | \textbf{s}_{N-1}, \textbf{A}) + \ln p(x_N | \textbf{s}_N, \textbf{λ})>_{q(\textbf{λ})q(\textbf{π})q(\textbf{A})q(\textbf{S}_{＼N})} + const \\
&=& <\ln p(\textbf{s}_N | \textbf{s}_{N-1}, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{N-1})} + <\ln p(x_N | \textbf{s}_N, \textbf{λ})>_{q(\textbf{λ})} + const \\
\end{eqnarray*}
これでひとまず、以下の分解近似の各項についてざっくり計算できた。
$$
p(\textbf{S}, \textbf{λ}, \textbf{A}, \textbf{π} | X) ≈ { \prod_{n=1}^{N} q(\textbf{s}_n) } q(\textbf{λ}, \textbf{A}, \textbf{π})
$$
さきほど求めた $\ln q(\textbf{λ}, \textbf{π}, \textbf{A})$ の近似式は以下である。
\ln q(\textbf{λ}, \textbf{π}, \textbf{A})
= \ln p(\textbf{λ}) + \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{S} | \textbf{π}, \textbf{A})>_{q(S)} + \sum_{n=1}^{N} <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{s})_n} + const
ここで、$\textbf{λ}$ に関する項が $\textbf{π},\textbf{A}$ の項と分離されていることがわかる。
よって、
\ln q(\textbf{λ}, \textbf{π}, \textbf{A})
= \ln q(\textbf{λ}) + \ln q(\textbf{π}, \textbf{A})
と書ける。ただし、
\ln q(\textbf{λ}) = \ln p(\textbf{λ}) + \sum_{n=1}^{N} <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{s})_n} + const
\ln q(\textbf{π}, \textbf{A})
= \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{S} | \textbf{π}, \textbf{A})>_{q(S)} + const
である。
$q(λ)$ を具体的に計算してみよう。
\begin{eqnarray*}
\ln q(\textbf{λ})
&=& \ln p(\textbf{λ}) + \sum_{n=1}^{N} <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{s})_n} + const \\
&=& \ln \prod_{k=1}^{K} Gam(\textbf{λ}_k | a,b) + \sum_{n=1}^{N}<\ln \prod_{k=1}^{K} Poi(x_n | λ_k)^{s_{n,k}}>_{q(\textbf{s})_n} + const \\
&=& \sum_{k=1}^{K} \ln Gam(\textbf{λ}_k | a,b) + \sum_{n=1}^{N}<\sum_{k=1}^{K} s_{n,k} \ln Poi(x_n | λ_k)>_{q(\textbf{s})_n} + const \\
&=& \sum_{k=1}^{K} \ln Gam(\textbf{λ}_k | a,b) + \sum_{n=1}^{N}\sum_{k=1}^{K} <s_{n,k}>_{q(\textbf{s})_n} \ln Poi(x_n | λ_k) + const \\
&=& \sum_{k=1}^{K} ((a-1) \ln λ_k - b λ_k + ln C_G(a,b)) + \sum_{n=1}^{N}\sum_{k=1}^{K} <s_{n,k}>_{q(\textbf{s})_n} (x_n \ln λ_k - \ln x_n! - λ_k) + const \\
&=& \sum_{k=1}^{K} \{ (a-1 + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})}) \ln λ_k - (b + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})})λ_k \} + const \\
\end{eqnarray*}
$\sum_{k=1}^{K}$ の表記から、各 $k$ について分解できることがわかる。つまり、
$$
q(\textbf{λ}) = \prod_{k=1}^{K} q(λ_k)
$$
である。
q(λ_k) = (a-1 + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})}) \ln λ_k - (b + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})})λ_k
はガンマ分布の対数表現であり、
q(λ_k) = Gam(λ_k | \hat{a_k}, \hat{b_k})
\hat{a_k} = a + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})}
\hat{b_k} = b + \sum_{n=1}^{N} <s_{n,k}>_{q(\textbf{S})}
であることがわかる。
実は $ q(\textbf{π}, \textbf{A}) $ も分解できる。
\begin{eqnarray*}
\ln q(\textbf{π}, \textbf{A})
&=& \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{S} | \textbf{π}, \textbf{A})>_{q(\textbf{S})} + const \\
&=& \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{s}_1 | \textbf{π}) \prod_{n=2}^{N} p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{S})} + const \\
&=& \ln p(\textbf{π}) + \ln p(\textbf{A}) + <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{s}_1)} + \sum_{n=2}^{N} <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{S})} + const \\
\end{eqnarray*}
よって、
\ln q(\textbf{π}, \textbf{A})
= \ln q(\textbf{π}) + \ln q(\textbf{A})
と書ける。ただし、
\ln q(\textbf{π})
= \ln p(\textbf{π}) + <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{s}_1)} + const
\ln q(\textbf{A})
= \ln p(\textbf{A}) + \sum_{n=2}^{N} <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{S})} + const
である。
よし、$q(\textbf{π})$ を具体的に計算してみよう。
\begin{eqnarray*}
\ln q(\textbf{π})
&=& \ln p(\textbf{π}) + <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{s}_1)} + const \\
&=& \ln Dir(\textbf{π}|\textbf{α}) + <\ln Cat(\textbf{s}_1 | \textbf{π})>_{q(\textbf{s}_1)} + const \\
&=& \ln C_D(\textbf{α}) \prod_{k=1}^{K} π_k^{α_k - 1} + <\ln \prod_{k=1}^{K} π_k^{s_1,k}>_{q(\textbf{s}_1)} + const \\
&=& \sum_{k=1}^{K} (α_k - 1) \ln π_k + \ln C_D(\textbf{α}) + < \sum_{k=1}^{K} s_{1,k} \ln π_k >_{q(\textbf{s}_1)} + const \\
&=& \sum_{k=1}^{K} (α_k - 1) \ln π_k + \ln C_D(\textbf{α}) + \sum_{k=1}^{K} <s_{1,k}>_{q(\textbf{s}_1)} \ln π_k + const \\
&=& \sum_{k=1}^{K} (α_k - 1 + <s_{1,k}>_{q(\textbf{s}_1)}) \ln π_k + const \\
\end{eqnarray*}
これはディリクレ分布の対数表現である。
q(\textbf{π}) = Dir(\textbf{π}|\hat{\textbf{α}})
ただし、
\hat{α}_k = <s_{1,k}>_{q(\textbf{s}_1)} + α_k
である。
次は $q(\textbf{A})$ を具体的に計算してみよう。
\begin{eqnarray*}
\ln q(\textbf{A})
&=& \ln p(\textbf{A}) + \sum_{n=2}^{N} <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{S})} + const \\
&=& \ln \prod_{i=1}^{K} Dir(\textbf{A}_{:,i} | \textbf{β}_{:,i}) + \sum_{n=2}^{N} <\ln \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j,i}^{\textbf{s}_{n,j} \textbf{s}_{n-1,i}}>_{q(\textbf{S})} + const \\
&=& \sum_{i=1}^{K} \ln Dir(\textbf{A}_{:,i} | \textbf{β}_{:,i}) + \sum_{n=2}^{N} < \sum_{i=1}^{K} \sum_{j=1}^{K} \textbf{s}_{n,j} \textbf{s}_{n-1,i} \ln A_{j,i}>_{q(\textbf{S})} + const \\
&=& \sum_{i=1}^{K} \ln C_D(\textbf{β}_{:,i}) \prod_{j=1}^{K} A_{j,i}^{β_{j,i} - 1} + \sum_{n=2}^{N} \sum_{i=1}^{K} \sum_{j=1}^{K} < \textbf{s}_{n,j} \textbf{s}_{n-1,i} \ln A_{j,i}>_{q(\textbf{S})} + const \\
&=& \sum_{i=1}^{K} \sum_{j=1}^{K} ((β_{j,i} - 1) \ln A_{j,i}) + \sum_{i=1}^{K} \sum_{j=1}^{K} \{ \sum_{n=2}^{N} < \textbf{s}_{n,j} \textbf{s}_{n-1,i}>_{q(\textbf{S})} \ln A_{j,i} \} + const \\
&=& \sum_{i=1}^{K} \{ \sum_{j=1}^{K} (β_{j,i} - 1 + \sum_{n=2}^{N} < \textbf{s}_{n,j} \textbf{s}_{n-1,i}>_{q(\textbf{S})}) \ln A_{j,i} \} + const \\
\end{eqnarray*}
$\sum_{i=1}^{K}$ から $K$ 個の近似分布の積に分解できることがわかる。
q(A_{:,i}) = Dir(\textbf{A}_{:,i}, \hat{\textbf{β}}_{:,i})
ただし、
\hat{β}_{j.i} = β_{j,i} + \sum_{n=2}^{N} < \textbf{s}_{n,j} \textbf{s}_{n-1,i}>_{q(\textbf{S})}
ちなみに、今のところは完全分解変分推論の仮定により時系列方向が独立なので、$< \textbf{s}\_{n,j} \textbf{s}\_{n-1,i}>\_{q(\textbf{S})} = < \textbf{s}\_{n,j}>\_{q(\textbf{S})} <\textbf{s}\_{n-1,i}>_{q(\textbf{S})}$ である。
最後に、$s_i$ を具体的に計算してみよう。上で計算した結果を再掲。
$i=1$ のとき
\ln q(\textbf{s}_1)
= <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{π})} + <\ln p(\textbf{s}_2 | \textbf{s}_1, \textbf{A})>_{q(\textbf{s}_2)q(\textbf{A})} + <\ln p(x_1 | \textbf{s}_1, \textbf{λ})>_{q(\textbf{λ})} + const
$i=n~(n=2,..,N-1)$ のとき
\ln q(\textbf{s}_n)
= <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n-1})} + <\ln p(\textbf{s}_{n+1} | \textbf{s}_n, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n+1})} + <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} + const
$i=N$ のとき
\ln q(\textbf{s}_N)
= <\ln p(\textbf{s}_N | \textbf{s}_{N-1}, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{N-1})} + <\ln p(x_N | \textbf{s}_N, \textbf{λ})>_{q(\textbf{λ})} + const
上記の3パターンに共通する $<\ln p(x\_i | \textbf{s}\_i, \textbf{λ})>_{q(\textbf{λ})}$ について、先に計算しておこう。
\begin{eqnarray*}
<\ln p(x_i | \textbf{s}_i, \textbf{λ})>_{q(\textbf{λ})}
&=& <\ln \prod_{k=1}^{K} Poi(x_i | λ_k)^{s_{i,k}} >_{q(\textbf{λ})} \\
&=& <\sum_{k=1}^{K} s_{i,k} \ln Poi(x_i | λ_k) >_{q(\textbf{λ})} \\
&=& \sum_{k=1}^{K} s_{i,k} <\ln Poi(x_i | λ_k) >_{q(λ_k)} \\
&=& \sum_{k=1}^{K} s_{i,k} <x_i \ln λ_k - \ln x_i ! - λ_k>_{q(λ_k)} \\
&=& \sum_{k=1}^{K} s_{i,k} (x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
\end{eqnarray*}
それでは、それぞれの $s_i$ の計算を進めよう。
$i=1$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_1)
&=& <\ln p(\textbf{s}_1 | \textbf{π})>_{q(\textbf{π})} + <\ln p(\textbf{s}_2 | \textbf{s}_1, \textbf{A})>_{q(\textbf{s}_2)q(\textbf{A})} + <\ln p(x_1 | \textbf{s}_1, \textbf{λ})>_{q(\textbf{λ})} + const \\
&=& <\ln Cat(\textbf{s}_1 | \textbf{π})>_{q(\textbf{π})} + <\ln \prod_{j=1}^{K} Cat(\textbf{s}_2 | \textbf{A}_{:,j})^{s_{1,j}}>_{q(\textbf{s}_2)q(\textbf{A})} + \sum_{k=1}^{K} s_{1,k} (x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& <\ln \prod_{k=1}^{K} π_k^{s_{1,k}}>_{q(\textbf{π})} + <\ln \prod_{j=1}^{K} \prod_{k=1}^{K} A_{j,k}^{s_{1,k} s_{2,j}}>_{q(\textbf{s}_2)q(\textbf{A})} + \sum_{k=1}^{K} s_{1,k} (x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& <\sum_{k=1}^{K} s_{1,k} \ln π_k>_{q(\textbf{π})} + <\ln \prod_{j=1}^{K} \prod_{k=1}^{K} A_{j,k}^{s_{1,j} s_{2,j}}>_{q(\textbf{s}_2)q(\textbf{A})} + \sum_{k=1}^{K} s_{1,k} (x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{k=1}^{K} s_{1,k} <\ln π_k>_{q(\textbf{π})} + \sum_{k=1}^{K} \sum_{j=1}^{K} s_{1,k} <s_{2,j}>_{q(\textbf{s}_2)} <\ln A_{j,k}>_{q(\textbf{A})} + \sum_{k=1}^{K} s_{1,k} (x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{k=1}^{K} s_{1,k} \{ <\ln π_k>_{q(\textbf{π})} + \sum_{j=1}^{K} <s_{2,j}>_{q(\textbf{s}_2)} <\ln A_{j,k}>_{q(\textbf{A})} + x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)} \} + const \\
\end{eqnarray*}
これはカテゴリ分布の対数表現である。
q(\textbf{s}_1) = Cat(\textbf{s}_1 | \textbf{η}_1)
ただし、
η_{1,k} \propto exp \{ <\ln π_k>_{q(\textbf{π})} + \sum_{j=1}^{K} <s_{2,j}>_{q(\textbf{s}_2)} <\ln A_{j,k}>_{q(\textbf{A})} + x_i <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)} \}
\sum_{k=1}^{K} η_{1,k} = 1
である。
$i=n~(n=2,..,N-1)$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_n)
&=& <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n-1})} + <\ln p(\textbf{s}_{n+1} | \textbf{s}_n, \textbf{A})>_{q(\textbf{A})q(\textbf{s}_{n+1})} + <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} + const \\
&=& <\ln \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j,i}^{s_{n,j} s_{n-1,i}}>_{q(\textbf{A})q(\textbf{s}_{n-1})} + <\ln \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j,i}^{s_{n+1,j} s_{n,i}}>_{q(\textbf{A})q(\textbf{s}_{n+1})} + \sum_{k=1}^{K} s_{n,k} (x_n <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{i=1}^{K} \sum_{j=1}^{K} s_{n,j} <s_{n-1,i}>_{q(\textbf{s}_{n-1})} <\ln A_{j,i}>_{q(\textbf{A})} + \sum_{i=1}^{K} \sum_{j=1}^{K} <s_{n+1,i}>_{q(\textbf{s}_{n+1})} s_{n,j} <\ln A_{j,i}>_{q(\textbf{A})} + \sum_{k=1}^{K} s_{n,k} (x_n <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{i=1}^{K} \sum_{k=1}^{K} s_{n,k} <s_{n-1,i}>_{q(\textbf{s}_{n-1})} <\ln A_{k,i}>_{q(\textbf{A})} + \sum_{k=1}^{K} \sum_{j=1}^{K} <s_{n+1,j}>_{q(\textbf{s}_{n+1})} s_{n,k} <\ln A_{j,k}>_{q(\textbf{A})} + \sum_{k=1}^{K} s_{n,k} (x_n <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{k=1}^{K} s_{n,k} \{ \sum_{i=1}^{K} <s_{n-1,i}>_{q(\textbf{s}_{n-1})} <\ln A_{k,i}>_{q(\textbf{A})} + \sum_{j=1}^{K} <s_{n+1,j}>_{q(\textbf{s}_{n+1})} <\ln A_{j,k}>_{q(\textbf{A})} + (x_n <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) \} + const \\
\end{eqnarray*}
これはカテゴリ分布の対数表現である。
q(\textbf{s}_n) = Cat(\textbf{s}_n | \textbf{η}_n)
ただし、
η_{n,k} \propto exp \{ \sum_{i=1}^{K} <s_{n-1,i}>_{q(\textbf{s}_{n-1})} <\ln A_{k,i}>_{q(\textbf{A})} + \sum_{j=1}^{K} <s_{n+1,j}>_{q(\textbf{s}_{n+1})} <\ln A_{j,k}>_{q(\textbf{A})} + (x_n <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) \}
\sum_{k=1}^{K} η_{n,k} = 1
である。
$i=N$ のとき
\begin{eqnarray*}
\ln q(\textbf{s}_N)
&=& <\ln \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j,i}^{s_{N,j} s_{N-1,i}}>_{q(\textbf{A})q(\textbf{s}_{N-1})} + <\ln p(x_N | \textbf{s}_N, \textbf{λ})>_{q(\textbf{λ})} + const \\
&=& \sum_{i=1}^{K} \sum_{k=1}^{K} s_{N,k} <s_{N-1,i}>_{q(\textbf{s}_{N-1})}<\ln A_{k,i}>_q(\textbf{A}) + \sum_{k=1}^{K} s_{N,k} (x_N <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)}) + const \\
&=& \sum_{k=1}^{K} s_{N,k} \{ \sum_{i=1}^{K} <s_{N-1,i}>_{q(\textbf{s}_{N-1})}<\ln A_{k,i}>_q(\textbf{A}) + x_N <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)} \}+ const \\
\end{eqnarray*}
これはカテゴリ分布の対数表現である。
q(\textbf{s}_N) = Cat(\textbf{s}_N | \textbf{η}_N)
ただし、
η_{n,k} \propto exp \{ \sum_{i=1}^{K} <s_{N-1,i}>_{q(\textbf{s}_{N-1})}<\ln A_{k,i}>_q(\textbf{A}) + x_N <\ln λ_k>_{q(λ_k)} - <λ_k>_{q(λ_k)} \}
\sum_{k=1}^{K} η_{N,k} = 1
である。
これまでの式に登場した期待値は解析的に得られるらしい。
私自身が解析的に導出できなかったものがいくつかあるので、これについてはまたの機会に。
<s_{n,i}> = η_{n,i}
<λ_i> = \frac{\hat{a}_i}{\hat{b}_i}
<\ln λ_i> = ψ(\hat{a}_i) - \ln \hat{b}_i
<\ln π_i> = ψ(\hat{a}_i) - ψ(\sum_{i'=1}^{K} \hat{α}_i)
<\ln A_{j,i}> = ψ(\hat{β}_{j,i}) - ψ(\sum_{j'=1}^{K} \hat{β}_{j',i})
構造化変分推論
完全分解変分推論では、時系列方向の独立性を仮定していた。ただ、時系列方向の相関をとらえたいので、分解を仮定しない推論方法がある。それが、構造化変分推論である。これこそ隠れマルコフモデルのマルコフ性を実現しているようでいいね。
このとき、事後分布の分解近似は以下のようになる。
\begin{eqnarray*}
p(\textbf{S}, \textbf{λ}, \textbf{π}, \textbf{A} | \textbf{X})
&≈& q(\textbf{S}) q(\textbf{λ}, \textbf{π}, \textbf{A}) \\
&=& q(\textbf{S}) q(\textbf{λ}) q(\textbf{π}) q(\textbf{A}) \\
\end{eqnarray*}
ただし、上述の議論から、$q(\textbf{λ}, \textbf{π}, \textbf{A}) = q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})$ とした。
はい変分法。
\begin{eqnarray*}
\ln q(\textbf{S})
&=& <\ln p(\textbf{S}, \textbf{λ}, \textbf{π}, \textbf{A}, \textbf{X})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} \\
&=& <\ln p(\textbf{π}) + \ln p(\textbf{A}) + \ln p(\textbf{λ}) + \ln p(\textbf{S} | \textbf{π}, \textbf{A}) + \ln p(\textbf{X} | \textbf{S}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} \\
&=& \sum_{n=1}^{N} <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} + <\ln p(\textbf{S} | \textbf{π}, \textbf{A})>_{q(\textbf{π}) q(\textbf{A})} + const \\
\end{eqnarray*}
構造化変分推論では、時系列方向の分解ができないので、解析的にパラメータの更新式を求めることは極めて困難。ただ、パラメータ($\textbf{π}$ とか $\textbf{λ}$ とか)の更新に必要な特定の数式だけに絞って求めることを目指す。
具体的には、
<s_{n,i}>
と
<s_{n-1, i} s_{n, j}>
である。
理想を言えば、$q(\textbf{S})$ の周辺化によって、下記の二つを求めたいところ。
q(\textbf{s}_n) = \sum_{S＼n} q(\textbf{S})
q(\textbf{s}_{n-1}, \textbf{s}_n) = \sum_{S＼ \{ n-1, n \}} q(\textbf{S})
ただ、周辺化のコストはすごい。各 $n$ について $\textbf{s}_n$ の次元が $K$ 個あるので、$K^N$ オーダーの計算が必要になる($i \in { 1,2,..,n-1,n+1,..,N } $ について $\textbf{s}_i$ を周辺化するので、正確な計算量は $K^{N-1}$)。これでは、データの数 $N$ が増えた時に対応できない。
そこで用いられる手法が、フォワードバックワードアルゴリズムである。
なにか大胆な近似をするわけではない。$\textbf{s}_n$ に着目しながら $q(\textbf{S})$ を変形するだけだ。
\begin{eqnarray*}
\ln q(\textbf{S})
&=& <\ln p(\textbf{S}, \textbf{λ}, \textbf{π}, \textbf{A}, \textbf{S})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} \\
&=& <\ln p(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1} | \textbf{π}, \textbf{A}, \textbf{λ}) + \ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A}) + \ln p(x_n | \textbf{s}_n, \textbf{λ}) + p(\textbf{S}_{n+1:N}, \textbf{X}_{n+1:N} | \textbf{s}_n, \textbf{A}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} + const \\
\end{eqnarray*}
ここで、例えば $ \textbf{S}_{1:n-1} $ は $ { \textbf{s}_1, \textbf{s}_2, .. , \textbf{s}_{n-1} } $ を意味し、$\textbf{X}_{n+1:N}$ は $ { x_{n+1}, x_{n+2}, .. , x_N } $ を意味している。
上式より、以下を得る。
\begin{eqnarray*}
q(\textbf{S})
&\propto& exp \{ <\ln p(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1} | \textbf{π}, \textbf{A}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} + <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{A})} + <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} + <p(\textbf{S}_{n+1:N}, \textbf{X}_{n+1:N} | \textbf{s}_n, \textbf{A}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{A})} \} \\
&=& \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n)
\end{eqnarray*}
ただし、
\tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) = exp \{ <\ln p(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1} | \textbf{π}, \textbf{A}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{π}) q(\textbf{A})} \}
\tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) = exp \{ <\ln p(\textbf{s}_n | \textbf{s}_{n-1}, \textbf{A})>_{q(\textbf{A})} \}
\tilde{p}(x_n | \textbf{s}_n) = exp \{ <\ln p(x_n | \textbf{s}_n, \textbf{λ})>_{q(\textbf{λ})} \}
\tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) = exp \{ <p(\textbf{S}_{n+1:N}, \textbf{X}_{n+1:N} | \textbf{s}_n, \textbf{A}, \textbf{λ})>_{q(\textbf{λ}) q(\textbf{A})} \}
である。
上式を用いて $ q(\textbf{s}_n) $ を計算してみよう。
\begin{eqnarray*}
q(\textbf{s}_n)
&=& \sum_{S＼n} q(\textbf{S}) \\
&\propto& \sum_{S＼n} \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \\
&=& \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{1:n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \sum_{S_{n+1:N}} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \\
&=& \textbf{f}(\textbf{s}_n)・\textbf{b}(\textbf{s}_n) \\
\end{eqnarray*}
ただし、
\textbf{f}(\textbf{s}_n) = \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{1:n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1})
\textbf{b}(\textbf{s}_n) = \sum_{S_{n+1:N}} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n)
ここで、$\sum_{\textbf{S}_{a:b}}$ の意味について言及しておく。これは、${ \textbf{s}_a, \textbf{s}_{a+1}, .., \textbf{s}_b }$ について周辺化するという意味である。
のちほど、$\sum_{\textbf{S}_{n}}$ のような表現が出てきて戸惑うかもしれないが、これの意味は「$\textbf{s}_{n}$について周辺化する」である。
また、上式の $\textbf{f}(\textbf{s}_n)$ と $\textbf{b}(\textbf{s}_n)$ がフォワードバックワードアルゴリズムの名前に対応しているのだろう。
上式変形で、一部ん？と私が思った箇所があったので、そこについて言及する。私が戸惑ったのは、下記式変形。
\sum_{S＼n} \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n)
= \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{1:n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \sum_{S_{n+1:N}} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n)
これに、(私目線で)わかりやすく途中式を追加してみる。
\begin{eqnarray*}
\sum_{S＼n} \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n)
&=& \sum_{S_{ \{ 1,..,n-1, ~ n+1,.., N \} }} \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \\
&=& \sum_{S_{ \{ 1:n-1 \} }} \{ \sum_{S_{ \{ n+1:N \} }} \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \} \\
&=& \sum_{S_{ \{ 1:n-1 \} }} \{ \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{ \{ n+1:N \} }} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \} \\
&=& \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{1:n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \sum_{S_{n+1:N}} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \\
\end{eqnarray*}
では、$\textbf{f}(\textbf{s}_n)$ と $\textbf{b}(\textbf{s}_n)$ の性質について、次の式で見ていこう。
$\textbf{f}(\textbf{s}_n)$ については以下。
\begin{eqnarray*}
\textbf{f}(\textbf{s}_n)
&=& \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{1:n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{S}_{1:n-1}, \textbf{X}_{1:n-1}) \\
&=& \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \tilde{p}(\textbf{s}_{n-1}, \textbf{x}_{n-1}) \sum_{S_{1:n-2}} \tilde{p}(\textbf{s}_{n-1} | \textbf{s}_{n-2}) \tilde{p}(\textbf{S}_{1:n-2}, \textbf{X}_{1:n-2}) \\
&=& \tilde{p}(x_n | \textbf{s}_n) \sum_{S_{n-1}} \tilde{p}(\textbf{s}_n | \textbf{s}_{n-1}) \textbf{f}(\textbf{s}_{n-1}) \\
\end{eqnarray*}
$\textbf{b}(\textbf{s}_n)$ については以下。
\begin{eqnarray*}
\textbf{b}(\textbf{s}_n)
&=& \sum_{S_{n+1:N}} \tilde{p}(\textbf{S}_{n+1:N}, \textbf{X}_{n+1, N} | \textbf{s}_n) \\
&=& \sum_{S_{n+1}} \tilde{p}(\textbf{s}_{n+1}, \textbf{x}_{n+1} | \textbf{s}_n) \sum_{S_{n+2:N}} \tilde{p}(\textbf{S}_{n+2:N}, \textbf{X}_{n+2, N} | \textbf{s}_n) \\
&=& \sum_{S_{n+1}} \tilde{p}(\textbf{s}_{n+1} | \textbf{s}_n) \tilde{p}(\textbf{x}_{n+1} | \textbf{s}_{n+1}) \sum_{S_{n+2:N}} \tilde{p}(\textbf{S}_{n+2:N}, \textbf{X}_{n+2, N} | \textbf{s}_n) \\
&=& \sum_{S_{n+1}} \tilde{p}(\textbf{s}_{n+1} | \textbf{s}_n) \tilde{p}(\textbf{x}_{n+1} | \textbf{s}_{n+1}) \textbf{b}(\textbf{s}_{n+1}) \\
\end{eqnarray*}
なるほど、$\textbf{f}(\textbf{s}_n)$ はインデックスが小さくなる方向に漸化式が書け、$\textbf{b}(\textbf{s}_n)$ はインデックスが大きくなる方向に漸化式が書けるのか。
$n=3, N=6$として具体的に書き下してみると、雰囲気がつかめてくる。
$\textbf{f}(\textbf{s}_3)$ については以下。
\begin{eqnarray*}
\textbf{f}(\textbf{s}_3)
&=& \tilde{p}(x_3 | \textbf{s}_3) \sum_{S_2} \tilde{p}(\textbf{s}_3 | \textbf{s}_2) \textbf{f}(\textbf{s}_2) \\
&=& \tilde{p}(x_3 | \textbf{s}_3) \sum_{S_2} \tilde{p}(\textbf{s}_3 | \textbf{s}_2) \tilde{p}(x_2 | \textbf{s}_2) \sum_{S_1} \tilde{p}(\textbf{s}_2 | \textbf{s}_1) \textbf{f}(\textbf{s}_1) \\
\end{eqnarray*}
$\textbf{b}(\textbf{s}_3)$ については以下。
\begin{eqnarray*}
\textbf{b}(\textbf{s}_3)
&=& \sum_{S_4} \tilde{p}(\textbf{s}_4 | \textbf{s}_3) \tilde{p}(\textbf{x}_4 | \textbf{s}_4) \textbf{b}(\textbf{s}_{4}) \\
&=& \sum_{S_4} \tilde{p}(\textbf{s}_4 | \textbf{s}_3) \tilde{p}(\textbf{x}_4 | \textbf{s}_4) \textbf{b}\sum_{S_5} \tilde{p}(\textbf{s}_5 | \textbf{s}_4) \tilde{p}(\textbf{x}_5 | \textbf{s}_5) \textbf{b}(\textbf{s}_{5}) \\
&=& \sum_{S_4} \tilde{p}(\textbf{s}_4 | \textbf{s}_3) \tilde{p}(\textbf{x}_4 | \textbf{s}_4) \textbf{b}\sum_{S_5} \tilde{p}(\textbf{s}_5 | \textbf{s}_4) \tilde{p}(\textbf{x}_5 | \textbf{s}_5) \textbf{b} \sum_{S_6} \tilde{p}(\textbf{s}_6 | \textbf{s}_5) \tilde{p}(\textbf{x}_6 | \textbf{s}_6) \textbf{b}(\textbf{s}_{6}) \\
\end{eqnarray*}
再帰的に行えるから実装はあまり難しくなさそうなうえ、$\textbf{s}_i$ と $\textbf{s}_j$ $(i ≠ j)$ の周辺化を別々に計算できるので、計算量が抑えられるメリットを享受できる。