Learning Problem: HMM Training
- Learning Problem:
, D是训练数据
- HMM的训练过程需要使用Forward/Backward算法用于求隐状态的期望
- HMM的训练过程也需要使用EM算法，根据F/B算法得到的隐状态期望以及D求关于
的最大似然估计
极大似然估计示例
假设隐状态已知，取值分别是{H， G}，且观测变量已知，取值分别是{1, 2, 3}, 共有四条数据，如下：
如何通过上面的示例数据计算最大似然估计呢？？？
1.估计π（隐状态的初始概率分布）
思路：根据上述示例先初始化π，只需要通过统计每个隐状态在序列中出现在首位的次数/总的序列个数，所以例如此处，可以计算得到：
2.估计A（转移概率矩阵）
3.估计B（发射概率矩阵）
极大似然估计总结：上述的计算过程，仅当已知隐藏状态序列时，以上最大似然估计才有效。 但是，事实并非如此，我们并不知道隐藏状态。 因此，我们需要找到另一种方法来估算过A，B。
EM算法（原理）
那么如何使用EM算法来对A， B， π估计呢？
- 初始化
.初始化所有概率相等或者随机初始化
- 基于1中计算每个过渡/发射使用频率的期望值。用计算结果来估计隐变量
，这里面用到了forward/backward算法
- 基于2中计算的隐变量重新估计概率
- 重复2,3直至收敛
概率方法
1.推导
idea：如果对于任意一个序列，我们知道每个t时刻，从状态i到状态j的概率，我们对该序列所有时刻的统计加和即可以得到从状态i到状态j的概率
即从数学表述，给定,在t时刻的状态为i，t+1时刻状态为j的概率, 其中
表示长度为T的观测序列
根据条件概率公式，有如下结论
,如果把
看成整体，可以推导
，进一步得到
所以对上式可以进一步改写成
表达方式一：
当然也可以通过另一种表达方式书写上面的式子,上式中即下式中的
, 上式中
即下式中的
,其中下式中忽略了
表达方式二：
也就是说不管是表达式1还是表达式2，我们都需要分别计算分子部分和分母部分，
- 1）分子部分是可以通过forward/backward算法进行计算的，即forward的
与
，因为是两个时刻所以还需要计算两个时刻的转移概率，以及
时刻的发射概率，那么最终得到的分子部分计算公式如下：
- 2）分母部分
是给定观测序列
以及模型参数
我们定义, 那么可以得到
注意：上述的只是对于某个观测序列而言的某一个时刻，我们需要对该观测序列所有的时刻都执行上面的操作，然后求和, 但是求和后需要进行归一化使其称为概率，我们知道了分子是从状态i到状态j,那么分母可以是状态i确定的情况下，下一时刻出现的不同的状态j
进一步推导可得：
- 3）分母的另一种解释
假设给定某个序列的概率为
注：下面的公式中用到了条件概率和联合概率公式的转换，注意区分
如果使用上述的表述重新表达,可以得到
通过公式（1）和公式（2）对比可以得到
也就是说两种角度都可以求得我们最终的目标
2.推导
给定隐状态时，生成
的概率可以用
表示,同时我们通过上面的分析已知，
在t时刻状态为j的概率为
那么要计算相当于求给定状态为j时的观测变量为
的概率，即如下表示
其中表示当t时刻为k时才为1，否则为0，即示性函数
EM算法（伪代码）
initialize A and B(random initialize or all equal)
-
while not convergence:
-
E-Step
-
M-Step
-
return A, B
代码实现
import pandas as pd
import numpy as np
def forward(V, a, b, initial_distribution):
alpha = np.zeros((V.shape[0], a.shape[0]))
alpha[0, :] = initial_distribution * b[:, V[0]]
for t in range(1, V.shape[0]):
for j in range(a.shape[0]):
# Matrix Computation Steps
# ((1x2) . (1x2)) * (1)
# (1) * (1)
alpha[t, j] = alpha[t - 1].dot(a[:, j]) * b[j, V[t]]
return alpha
def backward(V, a, b):
beta = np.zeros((V.shape[0], a.shape[0]))
# setting beta(T) = 1
beta[V.shape[0] - 1] = np.ones((a.shape[0]))
# Loop in backward way from T-1 to
# Due to python indexing the actual loop will be T-2 to 0
for t in range(V.shape[0] - 2, -1, -1):
for j in range(a.shape[0]):
beta[t, j] = (beta[t + 1] * b[:, V[t + 1]]).dot(a[j, :])
return beta
def em_learning(V, a, b, initial_distribution, n_iter=100):
M = a.shape[0]
T = len(V)
for n in range(n_iter):
alpha = forward(V, a, b, initial_distribution)
beta = backward(V, a, b)
xi = np.zeros((M, M, T - 1))
for t in range(T - 1):
denominator = np.dot(np.dot(alpha[t, :].T, a) * b[:, V[t + 1]].T, beta[t + 1, :])
for i in range(M):
numerator = alpha[t, i] * a[i, :] * b[:, V[t + 1]].T * beta[t + 1, :].T
xi[i, :, t] = numerator / denominator
gamma = np.sum(xi, axis=1)
a = np.sum(xi, 2) / np.sum(gamma, axis=1).reshape((-1, 1))
# Add additional T'th element in gamma
gamma = np.hstack((gamma, np.sum(xi[:, :, T - 2], axis=0).reshape((-1, 1))))
K = b.shape[1]
denominator = np.sum(gamma, axis=1)
for l in range(K):
b[:, l] = np.sum(gamma[:, V == l], axis=1)
b = np.divide(b, denominator.reshape((-1, 1)))
return {"a":a, "b":b}
data = pd.read_csv('data/data_python.csv')
V = data['Visible'].values
# Transition Probabilities
a = np.ones((2, 2))
A = a / np.sum(a, axis=1)
# Emission Probabilities
b = np.array(((1, 3, 5), (2, 4, 6)))
B = b / np.sum(b, axis=1).reshape((-1, 1))
# Equal Probabilities for the initial distribution
π = np.array((0.5, 0.5))
print(em_learning(V, A, B, π, n_iter=100))
{'a': array([[0.53816345, 0.46183655],
[0.48664443, 0.51335557]]), 'b': array([[0.16277513, 0.26258073, 0.57464414],
[0.2514996 , 0.27780971, 0.47069069]])}
