第4回 2019/5/13

講師：浅井潔

※ 今回の式にふってある式番号は、PRMLの対応する式番号です。

Definition of HMM

Markov Process（マルコフ連鎖）

系列データ $\mathbf{X} = \mathbf{x_1,\ldots,\mathbf{x}_N}$ があった時、未来の予測値が直近 $k$ 個の観測値以外の過去の観測値に対して独立であるという仮定をもち、条件付き分布 $p(\mathbf{x}_n|\mathbf{x}_{n-k},\ldots,\mathbf{x}_n-1)$ が $n$ によらずみな同一である時、この確率過程を定常 $k$ 次マルコフ連鎖と呼ぶ。

Markov Model（マルコフモデル）

マルコフモデルは、マルコフ連鎖の系列データ $\mathbf{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_N\}$ そのものを出力列としてみなす。（隠れマルコフモデルとの対応でこういう書き方をするが、当たり前）

例えば定常一次マルコフモデルの同時確率は以下で表される。

$$ p\left(\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}\right)=p\left(\mathbf{x}_{1}\right) \prod_{n=2}^{N} p\left(\mathbf{x}_{n} | \mathbf{x}_{n-1}\right)\qquad (13.2) $$

Hidden Markov Model（隠れマルコフモデル）

隠れマルコフモデルでは、マルコフ連鎖の系列データ $\mathbf{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_N\}$ を状態として考え、その状態からの出力記号列 $\mathbf{Z}=\{\mathbf{z}_1,\ldots,\mathbf{z}_N\}$ だけが観察できると考える。

このモデルは、状態（潜在変数）を 「$K$ 個の離散的な多項変数であり、$1$ 次のマルコフ性を有する」と仮定するのが一般的で、このモデルは以下のように表される。

$$ p\left(\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}, \mathbf{z}_{1}, \ldots, \mathbf{z}_{N}\right)=p\left(\mathbf{z}_{1}\right)\left[\prod_{n=2}^{N} p\left(\mathbf{z}_{n} | \mathbf{z}_{n-1}\right)\right] \prod_{n=1}^{N} p\left(\mathbf{x}_{n} | \mathbf{z}_{n}\right)\qquad (13.6) $$

この時、潜在変数は $K$ 次元の二値変数（one-of-K符号化法）なので、この条件付き分布は遷移確率(transition probability)を要素にもつ数表 $A$ に対応する。

つまり、遷移確率は $A_{jk}\equiv p(z_{nk}=1|z_{n-1,j}=1)$ で定義される。なお、確率なのでもちろん $0\ll A_{jk}\ll 1$ と $\sum_kA_{jk}=1$ を満たす。したがって、行列 $A$ は $K(K-1)$ 個の独立なパラメータを持つことになる。

この結果、条件付き分布は以下の形で明示的にかける。

$$ p\left(\mathbf{z}_{n} | \mathbf{z}_{n-1}, \mathbf{A}\right)=\prod_{k=1}^{K} \prod_{j=1}^{K} A_{j k}^{z_{n-1, j} z_{n k}}\qquad (13.7) $$

最初の潜在ノード $\mathbf{z}_1$ は、親ノードを持たないという点で特別である。このノードは、要素 $\pi_k\equiv p(z_{1k}=1)$ を持つ確率のベクトル $\boldsymbol{\pi}$ で表される周辺分布 $p(\boldsymbol{\pi})$ を持つ。$\left(\sum_k\pi_k=1\right)$

$$ p\left(\mathbf{z}_{1} | \boldsymbol{\pi}\right)=\prod_{k=1}^{K} \pi_{k}^{z_{1 k}}\qquad (13.8) $$

最後に、観測変数の条件付き確率分布 $p(\mathbf{x}_n|\mathbf{z}_n,\boldsymbol{\phi})$ を定義する。ここで、$\boldsymbol{\phi}$ は分布を支配するパラメータの集合である。

例えば $\mathbf{x}$ が連続変数の場合は出力各位率は以下で与えられる。$\left(\boldsymbol{\phi} = \{\boldsymbol{\mu}, \boldsymbol{\Sigma}\} \right)$

$$ p(\mathbf{x}_n | \mathbf{z}_n)=\prod_{k=1}^{K} \mathcal{N}\left(\mathbf{x}_n | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)^{z_{nk}}\qquad (9.11) $$

一方で $\mathbf{x}$ が離散的な場合は、条件付き確率表で与えられる。$\mathbf{x}_n$ は観測されるので、$\boldsymbol{\phi}$ の値が与えられたとき、分布 $p(\mathbf{x}_n|\mathbf{z}_n,\boldsymbol{\phi})$ は二値のベクトル $\mathbf{z}_n$ の $K$ 個の可能な状態に対応した $K$ 個の要素をもつベクトル（ベクトルのサイズは$\mathbf{x}$ の種類数）からなる。

以上をまとめると、同時確率は以下のように表される。

$$ \begin{aligned} p(\mathbf{X}, \mathbf{Z} | \boldsymbol{\theta}) &=p\left(\mathbf{z}_{1} | \boldsymbol{\pi}\right)\left[\prod_{n=2}^{N} p\left(\mathbf{z}_{n} | \mathbf{z}_{n-1}, \mathbf{A}\right)\right] \prod_{m=1}^{N} p\left(\mathbf{x}_{m} | \mathbf{z}_{m}, \boldsymbol{\phi}\right) & (13.10)\\ &=\prod_{k=1}^{K} \pi_{k}^{z_{1 k}} \left[\prod_{n=2}^{N} \prod_{i=1}^{K} \prod_{j=1}^{K} A_{j i}^{z_{n-1, j} z_{n i}}\right] \prod_{m=1}^{N}p\left(\mathbf{x}_{m} | \mathbf{z}_{m}, \boldsymbol{\phi}\right) \end{aligned} $$

小休止：Regular Grammar to HMM

正規文法(RG)の生成規則の全てに確率を付与することにより、確率正規文法(Stochastic Regular Grammar;SRG)となる。

確率正規文法は、非終端記号を状態、開始非終端記号を初期状態、終端記号を出力記号とみなすことにより、マルコフ連鎖から導入した隠れマルコフモデルと類似の確率モデルとなる。

ただし、確率正規文法では、各状態からではなく、各状態遷移から記号が出力されるという違いがある。

DP of HMM

Viterbi Algorithm

Viterbi アルゴリズムは、モデルパラメータ $\{\boldsymbol{\pi},\mathbf{A},\boldsymbol{\phi}\}$ がわかっているHMM及び配列 $\mathbf{X}$ が観測された時に、$\mathbf{X}$ を出力する最も確からしい状態系列 $\mathbf{Z}$ を求めるmax-sumアルゴリズムである。

$$ v_k(n) = p(\mathbf{x}_n|z_{nk}=1)\max_i\biggl\{v_i(n-1)p(z_{nk}=1|z_{n-1,i}=1)\biggr\}\\ v_k(1) = p(\mathbf{x}_1|z_{nk}=1) $$

ここで、$v_k(n)$ は、最後の状態 $\mathbf{z}_n$ が $z_{nk}=1$ であるという条件の下で $\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}$ を出力する最も確からしいパスの確率を表している。

動的計画法が終われば、あとは $\max_kv_k(N)$ の隠れ状態 $k$ からトレースバックすれば、求めたい状態系列 $\mathbf{Z}$ は求められる。

$K$ 状態からなるHMMと長さ $N$ の系列データ $\mathbf{X}$ が与えられた時、$\mathbf{X}$ の出力確率が最大となるパスは、Viterbi アルゴリズムにより、$O(NK^2)$ 時間で計算できる。

Forward Algorithm

Viterbi アルゴリズムでは、確率が最大となる1つのパス（上でいう赤いパス）のみを考えていたが、Forward アルゴリズムでは、$\mathbf{X}$ を出力する全ての確率を考えて、$\mathbf{X}$ が出力される確率を計算する。

この確率も、動的計画法を用いることで計算できる。ここで、$f_k(n)$ は、最後の状態 $\mathbf{z}_n$ が $z_{nk}=1$ であるという条件の下で $\{\mathbf{x}_1,\ldots,\mathbf{x}_n\}$ を出力する確率であるとする。

すると、この $f_k(n)$ は、次の再帰式に基づく動的計画法アルゴリズムにより計算することができる。

$$ f_k(n) = p(\mathbf{x}_n|z_{nk}=1)\sum_i^K\biggl\{f_i(n-1)p(z_{nk}=1|z_{n-1,i}=1)\biggr\}\\ f_k(1) = p(\mathbf{x}_1|z_{nk}=1) $$

※ Forward アルゴリズムとViterbi アルゴリズムの違いは、$\sum$ を計算しているか $\max$ を計算しているかだけである。

Backward Algorithm

Backward アルゴリズムは、出力系列 $\mathbf{X}$ が得られる確率を、Forward アルゴリズムとは逆向きの順序で計算する。

その目的は、HMMのパラメータ学習に用いられるBaum-Welch アルゴリズムなど、周辺確率の計算に用いることである。

次の再帰式を考える。

$$ b_k(n) = \sum_i^K\biggl\{p(\mathbf{x}_{n+1}|z_{n+1,i}=1)b_i(n+1) p(z_{n+1,i}=1|z_{nk}=1)\biggr\}\\ \forall_kb_k(N) = 1 $$

なお、Forward アルゴリズム、Backward アルゴリズムともに $O(NK^2)$ 時間で全ての $f_k(n),b_k(n)$ を計算できる。

EM of HMM

Marginalized Probabilities

States and Outputs

観測系列が $\mathbf{X} = \{\mathbf{x}_1,\ldots,\mathbf{x}_N\}$ であるとき、状態 $\mathbf{z}_n$ が $z_{nk}=1$ である確率は，

$$ \begin{aligned} p\left(z_{nk}=1, | \mathbf{X} \right) &=\frac{p(\mathbf{X}, z_{nk}=1)}{p(\mathbf{X})} \\ &=\frac{p\left(\mathbf{x}_1,\ldots,\mathbf{x}_n, z_{nk}=1 \right) p\left(\mathbf{x}_{n+1},\ldots,\mathbf{x}_N | z_{nk}=1 \right)}{p\left(\mathbf{X}\right)} \\ &=\frac{f_{k}(n) b_{k}(n)}{\sum_{i} f_{i}(T)} \end{aligned} $$

EM Algorithm

EMアルゴリズムは、観測できないデータ（隠れ変数）がある場合の最尤推定(Maximum likelihood estimation)のためのアルゴリズムで、次の式で定義される対数尤度(log likelihood)を最大化する $\boldsymbol{\theta}$ を計算することが目標となる。

$$\log p(\mathbf{X}|\boldsymbol{\theta}) = \log\sum_{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$$

この式の最適解を解析的に求めるのは困難であるので、反復を繰り返すことで対数尤度を単調に増加させ、その極大化を図る。

まず、上の式の左辺を以下のように変形する。（∵確率の乗法定理）

$$\log p(\mathbf{X}|\boldsymbol{\theta}) = \log p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) - \log p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})$$

EMアルゴリズムは、反復でパラメータを更新することによって対数尤度を増加させていく。そこで、一つ前のステップで求めたパラメータを $\boldsymbol{\theta}^{\mathrm{old}}$ と定義する。

すると、上の式の両辺に $p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})$ をかけて $\mathbf{Z}$ についての和をとると次の式を得る。

$$\log p(\mathbf{X}|\boldsymbol{\theta}) = \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})\log p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) -\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})\log p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})$$

右辺の第１項を $Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{\mathrm{old}})$ と定義する。つまり、以下のように表される。

$$ \begin{aligned} \log p(\mathbf{X}|\boldsymbol{\theta}) &= Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{\mathrm{old}}) -\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})\log p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})\\ \log p(\mathbf{X}|\boldsymbol{\theta}^{\mathrm{old}}) &= Q(\boldsymbol{\theta}^{\mathrm{old}}|\boldsymbol{\theta}^{\mathrm{old}}) -\sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})\log p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}}) \end{aligned} $$

したがって、以下の式が得られる。

$$ \begin{aligned} \log p(\mathbf{X}|\boldsymbol{\theta}) & - \log p(\mathbf{X}|\boldsymbol{\theta}^{\mathrm{old}})\\ & = Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{\mathrm{old}}) - Q(\boldsymbol{\theta}^{\mathrm{old}}|\boldsymbol{\theta}^{\mathrm{old}}) + \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})\log\frac{p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{\mathrm{old}})}{ p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})} \end{aligned} $$

ここで、最後の項はカルバックライブラー情報量(Kullback–Leibler divergence, 相対エントロピー)であり、常に非負である。なお、カルバックライブラー情報量の定義（離散）は以下で、 $P=Q$ で右辺の等号が成立する。

$$ \mathrm{KL}(P \| Q)=\sum_{i} P(i) \log \frac{P(i)}{Q(i)} \geq 0 $$

したがって、カルバックライブラー情報量が非負であることを踏まえると、以下の関係が成立する。

$$ \log p(\mathbf{X}|\boldsymbol{\theta}) - \log p(\mathbf{X}|\boldsymbol{\theta}^{\mathrm{old}}) \geq Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{\mathrm{old}}) - Q(\boldsymbol{\theta}^{\mathrm{old}}|\boldsymbol{\theta}^{\mathrm{old}}) $$

ゆえに、$\boldsymbol{\theta}^{\mathrm{new}} = \mathrm{argmax}_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{\mathrm{old}})$ とすることによって右辺が非負になり、尤度は増大するか、変化しないままとなる。

したがって、パラメータの更新を尤度関数が変化しなくなるまで行えば良い。

まとめ

初期パラメータ $\boldsymbol{\theta}^0$ を決定し、$t=0$ とする。
$Q(\boldsymbol{\theta}|\boldsymbol{\theta}^t) = \sum_{\mathbf{Z}}p(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^t)\log p(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$ を計算する。
$Q(\boldsymbol{\theta}|\boldsymbol{\theta}^t)$ を最大化する $\boldsymbol{\theta}^*$ を計算して $\boldsymbol{\theta}^{t+1} = \boldsymbol{\theta}^*$ とし、$t=t+1$ とする。
$Q(\boldsymbol{\theta}|\boldsymbol{\theta}^t)$ が増加しなくなるまでステップ $2,3$ を繰り返す。

Baum-Welch Algorithm

Baum-Welchアルゴリズムは、$f_k(n)$ と $b_k(n)$ を用いてHMMのパラメータの推定を行うEMアルゴリズムの一種である。

実際に配列解析に適用した例をここでやっているので、具体的な式変形はここで。

基本的には $Q$ 関数内で各パラメータが関与する部分が独立であるので、確率の制約条件（$[0,1]$ の範囲、足して $1$ になる）を制約条件として加えて、ラグランジュの未定乗数法を使うことになる。

ゲノム配列解析論Ⅰ 第4回

第4回 2019/5/13

Definition of HMM

Markov Process（マルコフ連鎖）

Markov Model（マルコフモデル）

Hidden Markov Model（隠れマルコフモデル）

小休止：Regular Grammar to HMM

DP of HMM

Viterbi Algorithm

Forward Algorithm

Backward Algorithm

EM of HMM

Marginalized Probabilities

States and Outputs

EM Algorithm

まとめ

Baum-Welch Algorithm

Table of Contents

Published

Last Updated

Category

Tags

Contact

Other contents