鮑姆-韋爾奇算法

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

歷史

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

介紹

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

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

記離散的隱含隨機變量為$X_{t}$，它總共有${\textstyle 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})$。

我們可以用${\textstyle \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$的概率。可以通過遞歸計算：

1. $\alpha _{i}(1)=\pi _{i}b_{i}(y_{1})$
2. $\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}$的概率。

1. $\beta _{i}(T)=1$
2. $\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。

例子

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

Transition State 1 State 2 State 1 0.5 0.5 State 2 0.3 0.7

Emission No Eggs Eggs State 1 0.3 0.7 State 2 0.8 0.2

Initial 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")，重新計算所有的轉移概率，得到下面的轉移矩陣：

Old Transition Matrix State 1 State 2 State 1 0.5 0.5 State 2 0.3 0.7

New Transition Matrix (Pseudo Probabilities) State 1 State 2 State 1 0.0598 0.0908 State 2 0.2179 0.9705

New Transition Matrix (After Normalization) 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

Old Emission Matrix No Eggs Eggs State 1 0.3 0.7 State 2 0.8 0.2

New Emission Matrix (Estimates) No Eggs Eggs State 1 0.0876 0.8769 State 2 1.0000 0.7385

New Emission Matrix (After Normalization) No Eggs Eggs State 1 0.0908 0.9092 State 2 0.5752 0.4248

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

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

代碼

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)))