鮑姆-韋爾奇演算法

在電氣工程、電腦科學、統計計算和生物資訊學中，鮑姆-韋爾奇演算法是用於尋找隱馬爾可夫模型未知參數的最大期望算法，它利用前向-後向演算法來計算E-Step的統計資訊。

## 歷史

鮑姆-韋爾奇演算法是以其發明者倫納德·埃紹·鮑姆和勞埃德·理察·韋爾奇的名字命名的。鮑姆-韋爾奇演算法和隱馬爾可夫模型在20世紀60年代末和70年代初由鮑姆和他的同事在國防分析研究所的一系列文章中首次描述。HMMs最初主要應用於語音處理領域。20世紀80年代，HMMs開始成為分析生物系統和資訊，特別是遺傳資訊的有用工具。此後，它們成為基因組序列機率建模的重要工具。

## 介紹

隱馬爾可夫模型描述了一組「隱含」變量和可觀測到的離散隨機變量的聯合機率。它依賴於假設：第$i$個隱藏變量只與第$i-1$個隱含變量相關，而與其他先前的隱藏變量無關，而當前觀測到的狀態僅依賴於當前的隱藏狀態。

鮑姆-韋爾奇算法利用最大期望算法，在給定一組觀測特徵向量的情況下，求出隱馬爾可夫模型參數的最大似然估計。

記離散的隱含隨機變量為$X_t$，它總共有$N$種狀態($X_t$有$N$個不同的值)。設$P(X_t|X_{t-1})$與時間無關，得到與時間無關的隨機變量轉移矩陣：

$A=\{a_{ij}\}=P(X_t=j|X_{t-1}=i)$初始的狀態(即$t=1$)分布由下式給出：

$$\pi_i=P(X_1=i)$$

記觀測到的變量為$Y_t$，總共有$K$種取值。同樣假設由隱含變量得到的可觀測變量與時間無關。在時間$t$，由隱含變量$X_t=j$得到的可觀察變量$Y_t=y_i$的機率是：

$$b_j(y_i)=P(Y_t=y_i|X_t=j)$$

由所有可能得$X_t$和$Y_t$的取值，我們可以得到$N\times K$的矩陣$B=\{b_j(y_i)\}$，其中$b_j$屬於所有可能得隱含狀態，$y_i$屬於所有的可觀測狀態。

給出可觀測序列：$Y=(Y_1=y_1, Y_2=y_2, \cdots, Y_T=y_T)$。

我們可以用$\theta(A, B, \pi)$描述隱馬爾科夫鏈，鮑姆-韋爾奇算法尋找$\theta^*=\arg \underset{\theta}{max}P(Y|\theta)$的局部極大值，也就是能夠使得觀測到的序列出現的機率最大的HMM的參數$\theta$。

## 演算法

初始化參數$\theta(A, B, \pi)$，可以隨機初始化，或者根據先驗知識初始化。

### 前向過程

記$\alpha_i(t)=P(Y_1=y_1, Y_2=y_2, \cdots, Y_t=y_t, X_t=i|\theta)$是參數$\theta$的條件下，觀測的序列是$y_1, y_2, \cdots, y_t$，時刻$t$的狀態是$i$的機率。可以通過遞歸計算：

- $$\alpha_i(1)=\pi_i b_i(y_1)$$
- $$\alpha_i(t+1) = b_i (y_{t+1})\sum_{j=1}^N\alpha_j(t)a_{ji}$$

### 後向過程

記$\beta_i(t) = P(Y_{t+1}=y_{t+1}, \cdots, Y_T=y_T|X_t=i, \theta)$是參數是$\theta$，在時刻$t$的狀態是$i$的條件下，餘下部分的觀測序列是$y_{t+1}, \cdots, y_{T}$的機率。

- $$\beta_i(T)=1$$
- $$\beta_i(t) = \sum_{j=1}^N \beta_j(t+1)a_{ij}b_j(y_{t+1})$$

### 更新

- 根據貝葉斯公式計算臨時變量。
- 在給定觀測序列$Y$和參數$\theta$的情況下，在時間$t$狀態是$i$的機率:$\gamma_i(t) = P(X_t=i|Y, \theta) = \frac{P(X_t=i, Y|\theta)}{P(Y|\theta)}=\frac{\alpha_i(t)\beta_i(t)}{\sum_{j=1}^N \alpha_j(t)\beta_j(t)}$
- 在給定觀測序列$Y$和參數$\theta$的情況下，在時間$t$狀態是$i$，在時間$t+1$狀態是$j$的機率:$\xi_{ij}(t) = P(X_t=i, X_{t+1}=j|Y, \theta) = \frac{P(X_t=i, X_{t+1}=j, Y|\theta)}{P(Y|\theta)}=\frac{\alpha_i(t)a_{ij}\beta_j(t+1)b_j(y_{t+1})}{\sum_{i=1}^N\sum_{j=1}^N\alpha_i(t) a_{ij} b_j(y_{t+1})\beta_j(t+1)}$
- $\gamma_i(t)$和$\xi_{ij}(t)$的分母一樣，表示給定參數$\theta$得到觀測序列$Y$的機率。
- $\pi_i^*=\gamma_i(1)$，在時間$1$狀態是$i$的機率
- $a_{ij}^* = \frac{\sum_{t=1}^{T-1}\xi_{ij}(t)}{\sum_{t=1}^{T-1}\gamma_i(t)}$，等於期望的從狀態$i$轉換到狀態$j$的數量除以從狀態$i$開始的轉換的總數。
- $b_i^*(v_k) = \frac{\sum_{t=1}^T 1_{y_t=v_k} \gamma_i(t)}{\sum_{t=1}^T\gamma_i(t)}$，其中$1_{y_t=v_k} = \begin{cases} 1\text{ if } y_t=v_k, \\ 0\text{ otherwise} \end{cases}$，$b_i^*(v_k)$是期望的從狀態$i$得到的觀察值等於$v_k$的數量除以從 狀態$i$開始的轉換的總數。
- 重複上面的步驟直到收斂。演算法可能過擬合，也不保證收斂到全局最大值。
- 其中計算$\gamma_i(t)$和$\xi_{ij}(t)$相當於最大期望算法的E-Step，而更新$\pi_i^* \alpha_{ij}^*, b_i^*(v_k)$的過程相當於最大期望算法的M-Step。

## 例子

假設我們有一隻會下蛋的雞，每天中午我們都會去拾取雞蛋。而雞是否下蛋依賴於一些未知的隱含狀態，這裡我們簡單的假設只有兩種隱含狀態會決定它是否下蛋。我們不知道這些隱含狀態的初始值，不知道他們之間的轉換機率，也不知道在每種狀態下雞會下蛋的機率。我們隨機初始化他們來開始猜測。

State 1｜State 2 State 1｜0.5｜0.5 State 2｜0.3｜0.7｜No Eggs｜Eggs State 1｜0.3｜0.7 State 2｜0.8｜0.2｜State 1｜0.2 State 2｜0.8

假設我們得到的觀測序列是(E=eggs, N=no eggs): N, N, N, N, N, E, E, N, N, N。

這樣我們同時也得到了觀測狀態的轉移：NN, NN, NN, NN, NE, EE, EN, NN, NN。

通過上面的資訊來重新估計狀態轉移矩陣。

Observed sequence｜Probability of sequence and state is S_1 then S_2｜Highest Probability of observing that sequence

NN｜0.024｜0.3584｜S_2,S_2

NN｜0.024｜0.3584｜S_2,S_2

NN｜0.024｜0.3584｜S_2,S_2

NN｜0.024｜0.3584｜S_2,S_2

NE｜0.006｜0.1344｜S_2,S_1

EE｜0.014｜0.0490｜S_1,S_1

EN｜0.056｜0.0896｜S_2,S_2

NN｜0.024｜0.3584｜S_2,S_2

NN｜0.024｜0.3584｜S_2,S_2

Total｜0.22｜2.4234

重新估計$S_1$到$S_2$的轉移機率為$\frac{0.22}{2.4234} = 0.0908$(下表中的"Pseudo probabilities")，重新計算所有的轉移機率，得到下面的轉移矩陣：

State 1｜State 2 State 1｜0.5｜0.5 State 2｜0.3｜0.7｜State 1｜State 2 State 1｜0.0598｜0.0908 State 2｜0.2179｜0.9705｜State 1｜State 2 State 1｜0.3973｜0.6027 State 2｜0.1833｜0.8167

接下來重新估計Emission Matrix:

Observed Sequence｜Highest probability of observing that sequence if E is assumed to come from S_1｜Highest Probability of observing that sequence

NE｜0.1344｜S_2,S_1｜0.1344｜S_2,S_1

EE｜0.0490｜S_1,S_1｜0.0490｜S_1,S_1

EN｜0.0560｜S_1,S_2｜0.0896｜S_2,S_2

Total｜0.2394｜0.2730

重新估計從隱含狀態$S_1$得到觀察結果E的機率是$\frac{0.2394}{0.2730} = 0.8769$，得到新的Emission Matrix

No Eggs｜Eggs State 1｜0.3｜0.7 State 2｜0.8｜0.2｜No Eggs｜Eggs State 1｜0.0876｜0.8769 State 2｜1.0000｜0.7385｜No Eggs｜Eggs State 1｜0.0908｜0.9092 State 2｜0.5752｜0.4248

為了估計初始狀態的機率，我們分別假設序列的開始狀態是$S_1$和$S_2$，然後求出最大的機率，再歸一化之後更新初始狀態的機率。

一直重複上面的步驟，直到收斂。

## 代碼

```python3
from hmmlearn import hmm
import numpy as np

X = np.array([1, 1, 1, 1, 1, 0, 0, 1, 1, 1]).reshape(-1, 1)
model = hmm.GaussianHMM(n_components=2, covariance_type='full')
model.fit(X)

model.monitor_.history

# pi
model.startprob_

# state transform matrix
model.transmat_

# emission_matrix
np.power(np.e, model._compute_log_likelihood(np.unique(X).reshape(-1, 1)))
```