隐马尔可夫模型的学习,根据训练数据是包括观测序列和状态序列还是只有观测序列,可以分别由监督学习与非监督学习实现。由于监督学习需要使用训练数据,而人工标注训练数据往往代价很高,有时就会利用非监督学习的方法,即Baum-Welch算法(也就是EM算法)。在介绍学习算法之前,先介绍一些概率和期望值的计算。这些计算会成为Baum-Welch算法公式的基础。
一些概率和期望值的计算
利用前向概率和后向概率,可以得到关于单个状态和两个状态概率的计算公式。
1. 给定模型
和观测
,在时刻
处于状态
的概率。记为
先分解为分数形式
根据前向概率的定义可以做以下变换
后向概率的定义如下
将这两者相乘得到
以上结果从两者的定义上也很好理解。
对变量 在范围 上求和
将式 代入 可以得到
2. 给定模型 和观测 ,在时刻 处于状态 且在时刻 处于状态 的概率。记为
通过前向后向概率计算:
分子可以用前向后向概率表示
则 可以表示为
3. 将 和 对各个时刻求和,可以得到一些有用的期望值。
(1) 观测 下,状态 出现的期望值
将每一个时刻下,出现状态 的概率相加
(2) 观测 下,由状态 转移的期望值
能够从状态 转移的时刻是 ,比上一个求和公式少了时刻
(3) 观测 下,由状态 转移到状态 的期望值
Baum-Welch模型
参数估计公式
·推导的过程,尤其是拉格朗日对偶,我暂时还不十分理解,先直接给出训练方法,公式和代码。Baum-Welch算法(Baum-Welch algorithm),它是EM算法在隐马尔可夫模型学习过程中的具体实现,由Baum和Welch提出。
(1)初始化
对n=0,选取
,得到模型
(2)递推。对
公式右端按照观测 和模型 代入计算
(3)终止,得到模型
Baum-Welch算法的Python实现
def baum_welch_train(self, observations, criterion=0.05):
n_states = self.A.shape[0]
n_samples = len(observations)
done = False
while not done:
# alpha_t(i) = P(O_1 O_2 ... O_t, q_t = S_i | hmm)
# Initialize alpha
alpha = self._forward(observations)
# beta_t(i) = P(O_t+1 O_t+2 ... O_T | q_t = S_i , hmm)
# Initialize beta
beta = self._backward(observations)
xi = np.zeros((n_states,n_states,n_samples-1))
for t in range(n_samples-1):
denom = np.dot(np.dot(alpha[:,t].T, self.A) * self.B[:,observations[t+1]].T, beta[:,t+1])
for i in range(n_states):
numer = alpha[i,t] * self.A[i,:] * self.B[:,observations[t+1]].T * beta[:,t+1].T
xi[i,:,t] = numer / denom
# gamma_t(i) = P(q_t = S_i | O, hmm)
gamma = np.sum(xi,axis=1)
# Need final gamma element for new B
prod = (alpha[:,n_samples-1] * beta[:,n_samples-1]).reshape((-1,1))
gamma = np.hstack((gamma, prod / np.sum(prod))) #append one more to gamma!!!
newpi = gamma[:,0]
newA = np.sum(xi,2) / np.sum(gamma[:,:-1],axis=1).reshape((-1,1))
newB = np.copy(self.B)
num_levels = self.B.shape[1]
sumgamma = np.sum(gamma,axis=1)
for lev in range(num_levels):
mask = observations == lev
newB[:,lev] = np.sum(gamma[:,mask],axis=1) / sumgamma
if np.max(abs(self.pi - newpi)) < criterion and \
np.max(abs(self.A - newA)) < criterion and \
np.max(abs(self.B - newB)) < criterion:
done = 1
self.A[:],self.B[:],self.pi[:] = newA,newB,newpi