Baum-Welch
Pair HMMとは、二本の配列を受け取ってペアワイズアラインメントを出力するHMMで、「アラインメントの確率的な解釈が可能になる」（アラインメントのどの部分が信頼できるかを確率的に評価できる）という利点があります。
このPair HMMは、一般的なHMMと同様にBaum-Welchアルゴリズムによってパラメータを最尤推定していきます。
\(\xi\)
HMM
一般的なHMMでは、例えば遷移確率 \(A_{jk}\)（状態 \(j\) から状態 \(k\) に遷移する確率）を最尤推定する際は、
- 直感的には「状態 \(j\) から状態 \(k\) に遷移する回数」を、「状態 \(j\) から遷移する回数」で割れば最尤解が求まる。
- 実際の隠れ状態の遷移が確定していないので、「期待される状態 \(j\) から状態 \(k\) に遷移する回数」を、「期待される状態 \(j\) から遷移する回数」で割れば良い。
- \(\mathbf{z}_{n-1}\) で状態 \(j\) におり、\(\mathbf{z}_n\) で状態 \(k\) にいる確率 \(\xi\left(z_{n-1,j},z_{nk}\right)\) は、forward-backward algorithmの \(\alpha,\beta\) を用いれば以下のように求めることができる。
$$\xi_n\left(\mathbf{z}_{n-1}, \mathbf{z}_{n}\right)=\frac{\alpha\left(\mathbf{z}_{n-1}\right) p\left(\mathbf{x}_{n} | \mathbf{z}_{n}\right) p\left(\mathbf{z}_{n} | \mathbf{z}_{n-1}\right) \beta\left(\mathbf{z}_{n}\right)}{p(\mathbf{X})}\qquad (13.43)$$
- \(\xi\) を用いて、最尤推定を行う。
$$A_{j k}= \frac{\sum_{n=2}^{N} \xi\left(z_{n-1, j}, z_{n k}\right)}{\sum_{l=1}^{K} \sum_{n=2}^{N} \xi\left(z_{n-1, j}, z_{n l}\right)} \qquad (13.19)$$
という様に行いますが、Pair HMMは、挿入・欠失を考えるため、\(n\) で配列 \(\mathbf{X,Y}\) のどの要素（\(\mathbf{x}_u,\mathbf{y}_v\)）を考えているのか一意に定まらないという点で普通のHMMとは異なります。
Pair HMM
そこで、「\(n\) ではなく \(u,v\) に注目する」ことで、 \((13.43)\) の式を以下のように書き換えます。
$$\xi_{u,v}\left(i,j\right)=\frac{f_i\left(u,v\right)A_{ij}\phi_j\left(\diamond_1\right)b_j\left(\diamond_2,\diamond_3\right)}{\sum_{k=1}^Kf_k\left(u,v\right)b_k\left(u,v\right)}$$
ここで、分子に注目すると、各変数の意味が
| variable | meaning |
|---|---|
| \(f_i\left(u,v\right)\) | The probability of ending up in state \(i\) after aligning the two sequences \(\mathbf{X,Y}\) up to observation \(u\) and \(v\) respectively. |
| \(A_{ij}\) | the transition probability from state \(i\) to state \(j\). |
| \(\phi_j\left(\diamond_1\right)\) | the emission probability of emitting \(\diamond_1\) in state \(j\). |
| \(b_j\left(\diamond_2,\diamond_3\right)\) | the probability of being in state \(j\), given the sequences \(\mathbf{X,Y}\) are aligned from observation \(\diamond_2+1\) and \(\diamond_3+1\) to the end of the sequences respectively. |
であるので、\(\diamond\) は状態によって異なり、具体的に書き下すと以下のようになります。
| state\variable | \(\phi_j\left(\diamond_1\right)\) | \(b_j\left(\diamond_2,\diamond_3\right)\) |
|---|---|---|
| \(M\) | \(e\left(x_{u+1},y_{v+1}\right)\) | \(b_j\left(u+1,v+1\right)\) |
| \(X\) | \(e\left(x_{u+1}\right)\) | \(b_j\left(u+1,v\right)\) |
| \(Y\) | \(e\left(y_{v+1}\right)\) | \(b_j\left(u,v+1\right)\) |
\(\gamma\)
同様に \(\gamma\) も考えると、
$$
\gamma_n\left(\mathbf{z}_{n}\right)= \frac{\alpha\left(\mathbf{z}_{n}\right) \beta\left(\mathbf{z}_{n}\right)}{p(\mathbf{X})} \qquad(13.33)$$
から、
$$\gamma_i\left(u,v\right)=\frac{f_i\left(u,v\right)b_j\left(u,v\right)}{\sum_{k=1}^Kf_k\left(u,v\right)b_k\left(u,v\right)}$$
のように書き換えることになります。
Maximization step
ここまでで \(\gamma,\xi\) が求まったので、M stepで各パラメータを更新します。なお、以下で \(w\) は全ての記号の組を表します。
$$w\in\left\{(k_x,k_y)|k_x,k_y\in\left\{A,T,G,C,\epsilon
\right\}\right\}\setminus(\epsilon,\epsilon)$$
\(\pi_k\)
$$\pi^{\star}_i = \sum_w\gamma_i^w(0,0)$$
\(A_{ij}\)
$$A^{\star}_{ij} = \frac{\sum_w\sum_{u=0}^{\tilde{U}}\sum_{v=0}^{\tilde{V}}\xi_{u,v}^{w}\left(i,j\right)}{\sum_w\sum_{u=0}^{\tilde{U}}\sum_{v=0}^{\tilde{V}}\sum_{j}^{K}\xi_{u,v}^{w}\left(i,j\right)}$$
| final emitting state | \(\tilde{U}\) | \(\tilde{V}\) |
|---|---|---|
| \(M\) | \(U-1\) | \(V-1\) |
| \(X\) | \(U-1\) | \(V\) |
| \(Y\) | \(U\) | \(V-1\) |
\(\phi_{i}(k)\)
$$\phi_i\left(k\right) = \frac{\sum_w\sum_{u=0}^{U\ast1}\sum_{v=0}^{V\ast2}\gamma_i^w\left(u,v\right)}{\sum_w\sum_{u=0}^U\sum_{v=0}^V\gamma_i^w\left(u,v\right)}$$
- \(\ast1\)：\(x_u=k_x\) and state \(i\) equals the state \(M\) or \(X\). (in the state \(Y\), a gap is present in observation stream \(x\) therefore \(k_x\) is not present.)
- \(\ast2\)：\(y_v=k_y\) and state \(j\) equals the state \(M\) or \(Y\).
Reference
- Wahle, Johannes and Armin Buch. "Alignment and word comparison with Pair Hidden Markov Models." (2013). pp.22-31
- Martijn B. Wieling. "Comparison of Dutch Dialects" (2007) pp.36-50