鲍姆-韦尔奇算法

在电气工程、计算机科学、统计计算和生物信息学中，鲍姆-韦尔奇算法是用于寻找隐马尔可夫模型未知参数的最大期望算法，它利用前向-后向算法来计算E-Step的统计信息。

历史

鲍姆-韦尔奇算法是以其发明者伦纳德·埃绍·鲍姆和劳埃德·理查德·韦尔奇的名字命名的。鲍姆-韦尔奇算法和隐马尔可夫模型在20世纪60年代末和70年代初由鲍姆和他的同事在国防分析研究所的一系列文章中首次描述。HMMs最初主要应用于语音处理领域。20世纪80年代，HMMs开始成为分析生物系统和信息，特别是遗传信息的有用工具。此后，它们成为基因组序列概率建模的重要工具。

介绍

隐马尔可夫模型描述了一组“隐含”变量和可观测到的离散随机变量的联合概率。它依赖于假设：第${\displaystyle i}$个隐藏变量只与第${\displaystyle i-1}$个隐含变量相关，而与其他先前的隐藏变量无关，而当前观测到的状态仅依赖于当前的隐藏状态。

鲍姆-韦尔奇算法利用最大期望算法，在给定一组观测特征向量的情况下，求出隐马尔可夫模型参数的最大似然估计。

记离散的隐含随机变量为${\displaystyle X_{t}}$，它总共有${\textstyle N}$种状态(${\displaystyle X_{t}}$有${\displaystyle N}$个不同的值)。设${\displaystyle P(X_{t}|X_{t-1})}$与时间无关，得到与时间无关的随机变量转移矩阵：

${\displaystyle A=\{a_{ij}\}=P(X_{t}=j|X_{t-1}=i)}$
初始的状态(即${\displaystyle t=1}$)分布由下式给出：

${\displaystyle \pi _{i}=P(X_{1}=i)}$

记观测到的变量为${\displaystyle Y_{t}}$，总共有${\displaystyle K}$种取值。同样假设由隐含变量得到的可观测变量与时间无关。在时间${\displaystyle t}$，由隐含变量${\displaystyle X_{t}=j}$得到的可观察变量${\displaystyle Y_{t}=y_{i}}$的概率是：

${\displaystyle b_{j}(y_{i})=P(Y_{t}=y_{i}|X_{t}=j)}$

由所有可能得${\displaystyle X_{t}}$和${\displaystyle Y_{t}}$的取值，我们可以得到${\displaystyle N\times K}$的矩阵${\displaystyle B=\{b_{j}(y_{i})\}}$，其中${\displaystyle b_{j}}$属于所有可能得隐含状态，${\displaystyle y_{i}}$属于所有的可观测状态。

给出可观测序列：${\displaystyle Y=(Y_{1}=y_{1},Y_{2}=y_{2},\cdots ,Y_{T}=y_{T})}$。

我们可以用${\textstyle \theta (A,B,\pi )}$描述隐马尔科夫链，鲍姆-韦尔奇算法寻找${\displaystyle \theta ^{*}=\arg {\underset {\theta }{max}}P(Y|\theta )}$的局部极大值，也就是能够使得观测到的序列出现的概率最大的HMM的参数${\displaystyle \theta }$。

算法

初始化参数${\displaystyle \theta (A,B,\pi )}$，可以随机初始化，或者根据先验知识初始化。

前向过程

记${\displaystyle \alpha _{i}(t)=P(Y_{1}=y_{1},Y_{2}=y_{2},\cdots ,Y_{t}=y_{t},X_{t}=i|\theta )}$是参数${\displaystyle \theta }$的条件下，观测的序列是${\displaystyle y_{1},y_{2},\cdots ,y_{t}}$，时刻${\displaystyle t}$的状态是${\displaystyle i}$的概率。可以通过递归计算：

1. ${\displaystyle \alpha _{i}(1)=\pi _{i}b_{i}(y_{1})}$
2. ${\displaystyle \alpha _{i}(t+1)=b_{i}(y_{t+1})\sum _{j=1}^{N}\alpha _{j}(t)a_{ji}}$

后向过程

记${\displaystyle \beta _{i}(t)=P(Y_{t+1}=y_{t+1},\cdots ,Y_{T}=y_{T}|X_{t}=i,\theta )}$是参数是${\displaystyle \theta }$，在时刻${\displaystyle t}$的状态是${\displaystyle i}$的条件下，余下部分的观测序列是${\displaystyle y_{t+1},\cdots ,y_{T}}$的概率。

1. ${\displaystyle \beta _{i}(T)=1}$
2. ${\displaystyle \beta _{i}(t)=\sum _{j=1}^{N}\beta _{j}(t+1)a_{ij}b_{j}(y_{t+1})}$

更新

• 根据贝叶斯公式计算临时变量。
  • 在给定观测序列${\displaystyle Y}$和参数${\displaystyle \theta }$的情况下，在时间${\displaystyle t}$状态是${\displaystyle i}$的概率:${\displaystyle \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)}}}$
  • 在给定观测序列${\displaystyle Y}$和参数${\displaystyle \theta }$的情况下，在时间${\displaystyle t}$状态是${\displaystyle i}$，在时间${\displaystyle t+1}$状态是${\displaystyle j}$的概率:${\displaystyle \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)}}}$
  • ${\displaystyle \gamma _{i}(t)}$和${\displaystyle \xi _{ij}(t)}$的分母一样，表示给定参数${\displaystyle \theta }$得到观测序列${\displaystyle Y}$的概率。

• 然后更新参数：
  • ${\displaystyle \pi _{i}^{*}=\gamma _{i}(1)}$，在时间${\displaystyle 1}$状态是${\displaystyle i}$的概率
  • ${\displaystyle a_{ij}^{*}={\frac {\sum _{t=1}^{T-1}\xi _{ij}(t)}{\sum _{t=1}^{T-1}\gamma _{i}(t)}}}$，等于期望的从状态${\displaystyle i}$转换到状态${\displaystyle j}$的数量除以从状态${\displaystyle i}$开始的转换的总数。
  • ${\displaystyle b_{i}^{*}(v_{k})={\frac {\sum _{t=1}^{T}1_{y_{t}=v_{k}}\gamma _{i}(t)}{\sum _{t=1}^{T}\gamma _{i}(t)}}}$，其中${\displaystyle 1_{y_{t}=v_{k}}={\begin{cases}1{\text{ if }}y_{t}=v_{k},\\0{\text{ otherwise}}\end{cases}}}$，${\displaystyle b_{i}^{*}(v_{k})}$是期望的从状态${\displaystyle i}$得到的观察值等于${\displaystyle v_{k}}$的数量除以从 状态${\displaystyle i}$开始的转换的总数。
• 重复上面的步骤直到收敛。算法可能过拟合，也不保证收敛到全局最大值。
• 其中计算${\displaystyle \gamma _{i}(t)}$和${\displaystyle \xi _{ij}(t)}$相当于最大期望算法的E-Step，而更新${\displaystyle \pi _{i}^{*}\alpha _{ij}^{*},b_{i}^{*}(v_{k})}$的过程相当于最大期望算法的M-Step。

例子

假设我们有一只会下蛋的鸡，每天中午我们都会去拾取鸡蛋。而鸡是否下蛋依赖于一些未知的隐含状态，这里我们简单的假设只有两种隐含状态会决定它是否下蛋。我们不知道这些隐含状态的初始值，不知道他们之间的转换概率，也不知道在每种状态下鸡会下蛋的概率。我们随机初始化他们来开始猜测。

假设我们得到的观测序列是(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 ${\displaystyle S_{1}}$ then ${\displaystyle S_{2}}$	Highest Probability of observing that sequence
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NE	0.006	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$
EE	0.014	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$
EN	0.056	0.0896	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
NN	0.024	0.3584	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
Total	0.22	2.4234

重新估计${\displaystyle S_{1}}$到${\displaystyle S_{2}}$的转移概率为${\displaystyle {\frac {0.22}{2.4234}}=0.0908}$(下表中的"Pseudo probabilities")，重新计算所有的转移概率，得到下面的转移矩阵：

接下来重新估计Emission Matrix:

Observed Sequence	Highest probability of observing that sequence if E is assumed to come from ${\displaystyle S_{1}}$		Highest Probability of observing that sequence
NE	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$	0.1344	${\displaystyle S_{2}}$,${\displaystyle S_{1}}$
EE	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$	0.0490	${\displaystyle S_{1}}$,${\displaystyle S_{1}}$
EN	0.0560	${\displaystyle S_{1}}$,${\displaystyle S_{2}}$	0.0896	${\displaystyle S_{2}}$,${\displaystyle S_{2}}$
Total	0.2394		0.2730

重新估计从隐含状态${\displaystyle S_{1}}$得到观察结果E的概率是${\displaystyle {\frac {0.2394}{0.2730}}=0.8769}$，得到新的Emission Matrix

为了估计初始状态的概率，我们分别假设序列的开始状态是${\displaystyle S_{1}}$和${\displaystyle 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)))