隐马尔可夫模型之Baum-Welch算法详解
前言
在上篇博文中，我们学习了隐马尔可夫模型的概率计算问题和预测问题，但正当要准备理解学习问题时，发现学习问题中需要EM算法的相关知识，因此，上一周转而学习了EM算法和极大似然估计，对隐藏变量的求解有了一些自己的理解，现在我们继续回过头来学习隐马尔可夫模型的学习问题。EM算法的相关介绍可参照博文 EM算法及其推广学习笔记。如果对隐马尔可夫模型还不胜了解的话，可参看博文隐马尔可夫学习笔记（一）。
学习问题
隐马尔可夫模型的学习，根据训练数据是包括观测序列和对应的状态序列还是只有观测序列，可以分别由监督学习与非监督学习实现。本节首先介绍监督学习算法，而后介绍非监督学习算法——Baum-Welch算法（也就是EM算法）。
监督学习问题
假设已给训练数据包含S个长度相同的观测序列和对应的状态序列{(O1,I1),(O2,I2),...,(OS,IS)}，那么可以利用极大似然估计方法来估计隐马尔可夫模型的参数，具体方法如下。
理论分析总算完毕了，简单总结一下前向后向算法，首先隐马尔可夫模型参数的估计问题是一个隐藏变量的极大似然估计，因此我们用到了EM算法来解决上述参数估计问题，从EM算法中，求得Q函数，从而能够对Q函数进行求偏导，得到极大似然函数的极值，求偏导算出了参数估计的公式，与先前λ^参数产生了关系，并进一步需要计算大量联合概率，而联合概率的计算巧妙的使用了节点图的各种性质，用中间变量降低了节点计算的复杂度，导出了对计算有帮助的定义，方便参数求解。
Code Time
Baum-Welch算法的Python实现
接着前文hmm.py和test.py文件继续添加。
1.在hmm.py中添加baum_welch_train算法
def baum_welch_train(self, observations, criterion=0.05):
n_states = self.A.shape[0]
# 观察序列的长度T
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_t(i)
alpha = self._forward(observations)
# beta_t(i) = P(o_t+1,o_t+2,...,o_T | q_t = s_i , hmm)
# Initialize beta
# 获得所有后向传播节点值 beta_t(i)
beta = self._backward(observations)
# 计算 xi_t(i,j) -> xi(i,j,t)
xi = np.zeros((n_states, n_states, n_samples - 1))
# 在每个时刻
for t in range(n_samples - 1):
# 计算P(O | hmm)
denom = sum(alpha[:, -1])
for i in range(n_states):
# numer[1,:] = 行向量，alpha[i,t]=实数，slef.A[i,:] = 行向量
# self.B[:,observations[t+1]].T = 行向量,beta[:,t+1].T = 行向量
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) 就是对j进行求和
gamma = np.sum(xi, axis=1)
# need final gamma elements for new B
prod = (alpha[:, n_samples - 1] * beta[:, n_samples - 1]).reshape((-1, 1))
# 合并T时刻的节点
gamma = np.hstack((gamma, prod / np.sum(prod)))
# 列向量
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, newpi2.在hmm.py中添加模拟序列生成函数
def simulate(self, T):
def draw_from(probs):
# np.random.multinomial 为多项式分布，1为实验次数，类似于投掷一枚骰子，丢出去是几，probs每个点数的概率，均为1/6
# 给定行向量的概率，投掷次数为1次，寻找投掷的点数
return np.where(np.random.multinomial(1, probs) == 1)[0][0]
observations = np.zeros(T, dtype=int)
states = np.zeros(T, dtype=int)
states[0] = draw_from(self.pi)
observations[0] = draw_from(self.B[states[0], :])
for t in range(1, T):
states[t] = draw_from(self.A[states[t - 1], :])
observations[t] = draw_from(self.B[states[t], :])
return observations, states回到海藻模型，我们可以用这样一串代码完成Baum-Welch算法的训练，并且评估其准确率。
3.在test.py中添加测试代码
# 参数估计
observations_data, states_data = h.simulate(100)
guess = hmm.HMM(array([[0.33, 0.33, 0.34],
[0.33, 0.33, 0.34],
[0.33, 0.33, 0.34]]),
array([[0.25, 0.25, 0.25, 0.25],
[0.25, 0.25, 0.25, 0.25],
[0.25, 0.25, 0.25, 0.25]]),
array([0.7, 0.15, 0.15])
)
guess.baum_welch_train(observations_data)
# 预测问题
states_out = guess.state_path(observations_data)[1]
p = 0.0
for s in states_data:
if next(states_out) == s: p += 1
print(p / len(states_data))经过多次测试，本算法的预测准确率在0.3~0.5。可见隐马尔可夫模型的参数估计的准确率还没有到令人满意的程度。
参考文献
- EM算法及其推广学习笔记
- 隐马尔可夫学习笔记（一）
- 李航. 统计学习方法[M]. 北京：清华大学出版社，2012