文章测试 HMM
1. 隐马尔可夫模型形式定义
所有可能的状态集合\(Q = \{q_1, q_2,...,q_N\}\)，所有可能的观测集合\(V = \{v_1, v_2,...,v_M\}\)
长度为\(T\)的状态序列\(I = (i_1, i_2,...,i_T)\)，\(I\)对应的观测序列\(O = (o_1, o_2,...,o_T)\)
状态转移矩阵$A = [a_{ij}]{N\times M} \(，\)a = P(i_{t+1}=q_j|i_t=q_i)\(，\)i=1,2,...,N\(，\)j=1,2,...,N\(
观测概率矩阵\)B = [b_j(k)]_{N\times M}\(，\)b_j(k) = P(o_t=v_k|i_t=q_j)\(，\)k=1,2...,M\(，\)j=1,2,...,N\(
初始状态概率向量\)\pi = (\pi_i)\(，\)\pi_i = P(i_1=q_i)\(，\)i=1,2,...,N\(
隐马尔可夫模型\)\lambda$可用三元式表示
\[\lambda = (\pi,A,B)
\]
2. 两个基本假设
- 齐次马尔可夫假设：假设隐藏的马尔可夫链在任意时刻\(t\)的状态只依赖于前一时刻的状态，与其他时刻状态、观测都无关，也与时刻\(t\)无关
\[P(i_t|i_1,o_1, i_2,o_2,...,i_{t-1},o_{t-1}) = P(i_t|i_{t-1})
\]
- 观测独立性假设：假设任意时刻的观测只与当前状态有关，与其他观测及状态无关
\[P(o_t|i_1,o_1, i_2,o_2,...,i_T,o_T) = P(o_t|i_t)
\]
3. 三个基本问题
- 概率计算问题
给定模型\(\lambda = (\pi, A, B)\)、观测序列\(O = (o_1, o_2,...,o_T)\)，计算模型\(\lambda\)下，\(O\)出现的概率\(P(O|\lambda)\)。 - 学习问题
已知观测序列\(O = (o_1, o_2,...,o_T)\)，用极大似然估计的方法估计模型\(\lambda = (\pi, A, B)\)的参数。 - 预测问题
已知模型\(\lambda = (\pi, A, B)\)和观测序列\(O = (o_1, o_2,...,o_T)\)，求最有可能的对应的状态序列\(I = (i_1, i_2,...,i_T)\)，\(P=(I|O)\)。
4. 概率计算问题
4.1 直接计算
计算所有可能的状态序列\(P(I|\lambda)\)，计算在状态序列为\(I\)时观测为\(O\)的概率\(P(O|I,\lambda)\)
\[P(O|\lambda) = \sum_IP(O|I,\lambda)P(I|\lambda)
\]
\(I\)和\(O\)的长度为\(T\)，\(I\)有\(N\)种取值，计算复杂度为\(O(N^TT)\)，实际不可行
4.2 前向算法
在\(t\)时刻，观测序列为\(o_1,o_2...,o_t\)下，状态为\(q_i\)的概率设为\(\alpha_t(i)=P(o_1,o_2,...,o_t,i_t=q_i|\lambda)\)
- 初值
\[\alpha_1(i) = \pi_ib_i(o_1) ,\quad i=1,2,...,N
\]
- 递推
\[\alpha_{t+1}(i) = \left (\sum_j^N \alpha_t(j)a_{ji} \right )b_i(o_{t+1})
\]
- 终止
\[P(O|\lambda) = \sum_i^N\alpha_{T}(i)
\]
每次计算\(N\)个状态情况，每个情况\(N\)次乘法，共\(T\)个时间步，计算复杂度为\(O(N^2T)\)
4.3 后向算法
在\(t\)时刻，当前状态\(i_t=q_i\)的情况下，后续的观测为\(o_{t+1},o_{t+2},...,o_T\)的概率\(\beta_t(i)= P(o_{t+1},o_{t+2},...,o_T|i_t=q_i,\lambda)\)
- 初值
\[\beta_T(i) = 1 ,\quad i=1,2,...,N
\]
- 递推
\[\beta_t(i) = \sum_j^Na_{ij}\beta_{t+1}(j)b_i(o_{t+1})
\]
- 终止
\[P(O|\lambda) = \sum_i^N\pi_i\beta_1(i)b_i(o_1)
\]
每次计算\(N\)个状态情况，每个情况\(2N\)次乘法，共\(T\)个时间步，计算复杂度为\(O(N^2T)\)
5. 学习问题
5.1 已知状态序列和对应的观测序列
给定训练集包含\(S\)个长度为\(T\)的状态序列和对应的观测序列\(\{(I_1, O_1),(I_2, O_2),...,(I_S, O_S)\}\)，那么可用极大似然法估计模型参数。
- 转移概率\(a_{ij}\)的估计
\[a_{ij} = \frac{1}{S}\sum_s^S \frac{\sum_i^{T-1}I(i_{s,i}=q_i,i_{s,i+1}=q_j)}{\sum_i^{T-1}I(i_{s,i}=q_i)}
\]
- 观测概率\(b_{ij}\)的估计
\[b_{ij} = \frac{1}{S}\sum_s^S \frac{\sum_i^{T}I(i_{s,i}=q_i,o_{s,i}=v_j)}{\sum_i^{T}I(i_{s,i}=q_i)}
\]
- 初始状态概率\(\pi_i\)的估计
\[\pi_i = \frac{1}{S}\sum_s^S I(i_{s,1} = q_i)
\]
5.2 已知观测序列
给定训练集包含\(S\)个长度为\(T\)的观测序列\(\{O_1,O_2,...,O_S\}\)。此时可将观测序列看作观测数据，将状态序列看作隐藏数据，该问题为含有隐藏变量的学习问题，可使用EM算法求解。
\[P(O|\lambda) = \sum_IP(O|I,\lambda)P(I|\lambda)
\]
- 选择\(\lambda\)的初值\(a_{ij}^{(0)}\)，\(b_{ij}^{(0)}\)，\(\pi_i^{(o)}\)
- E步
\[Q(\lambda|\overline \lambda ) = \sum_I \log P(O,I|\lambda)P(O,I|\overline \lambda)
\]
\[P(O,I|\lambda) = \pi_{i_1}b_{i_1}(o_1)\prod_{t=2}^{T}a_{i_{t-1}i_t}b_{i_t}(o_t)
\]
\[Q(\lambda|\overline \lambda) = \sum_I \log \pi_{i_1}P(O,I|\overline \lambda) +
\sum_I \left( \sum_{t=1}^T \log b_{i_t}(o_t)P(O,I|\overline \lambda) \right) + \\
\sum_I \left( \sum_{t=1}^{T-1} \log a_{i_ti_{t+1}} P(O,I|\overline \lambda)\right) \]
- M步
（1）对于第一项，\(\sum_{i=1}^N \pi_i = 1\)，引入拉格朗日乘子
\[\sum_I \log \pi_{i_1}P(O,I|\overline \lambda) + \gamma\left(\sum_{i=1}^N\pi_i-1\right)
\]
\[\sum_{i=1}^N \log \pi_i P(O,i_1=i|\overline \lambda) + \gamma\left(\sum_{i=1}^N\pi_i-1\right)
\]
对\(\pi_i\)求偏导令结果为0，得到
\[\frac{1}{\pi_i} P(O,i_1=i|\overline \lambda) + \gamma = 0
\]
\[P(O,i_1=i|\overline \lambda) + \gamma \pi_i =0
\]
进行求和
\[\sum_{i=1}^N\left( P(O,i_1=i|\overline \lambda) + \gamma \pi_i \right)=0
\]
\[P(O|\overline \lambda)=-\gamma
\]
最后得到
\[\pi_i = \frac{P(O,i_1=i|\overline \lambda)}{P(O|\overline \lambda)}
\]
（2）对于第二项
\[\sum_I \left( \sum_{t=1}^T \log b_{i_t}(o_t)P(O,I|\overline \lambda) \right)
\]
\[\sum_{i=1}^N \left( \sum_{t=1}^T \log b_i(o_t)P(O,i_t=i|\overline \lambda)\right)
\]
约束条件为
\[\sum_{k=1}^M b_i(k) = 1
\]
引入拉格朗日乘子
\[\sum_{i=1}^N \left( \sum_{t=1}^T \log b_i(o_t)P(O,i_t=i|\overline \lambda)\right) + \gamma \left( \sum_{k=1}^M b_i(k) -1\right)
\]
对\(b_i(k)\)求偏导，当\(o_t=v_k\)时\(b_i(k)\)对\(b_i(o_t)\)求导不为0，用\(I(o_t=v_k)\)表示。令结果为0，得到
\[\sum_{t=1}^T \frac{1}{b_i(k)}I(o_t=v_k)P(O,i_t=i|\overline \lambda)+\gamma=0
\]
\[\sum_{t=1}^T I(o_t=v_k)P(O,i_t=i|\overline \lambda)+\gamma b_i(k)=0
\]
进行求和
\[\sum_{k=1}^M \left( \sum_{t=1}^T I(o_t=v_k)P(O,i_t=i|\overline \lambda)+\gamma b_i(k)\right)=0
\]
对于一个时刻\(t\)输出\(o_t\)是确定的，在所有\(k=1,2,...,M\)的情况下，\(o_t=v_k\)的情况成立只有一个。也就是对\(k\)求和后，\(I(o_t=v_k)=1\)的情况有且仅有一个，得到
\[\sum_{t=1}^TP(O,i_t=i|\overline \lambda) + \gamma = 0
\]
\[\gamma = - \sum_{t=1}^TP(O,i_t=i|\overline \lambda)
\]
\[b_i(k) = \frac{ \sum\limits_{t=1}^T I(o_t=v_k)P(O,i_t=i|\overline \lambda)}{\sum\limits_{t=1}^TP(O,i_t=i|\overline \lambda)}
\]
（3）对于第三项
\[\sum_I \left( \sum_{t=1}^{T-1} \log a_{i_ti_{t+1}} P(O,I|\overline \lambda)\right)
\]
\[\sum_{i=1}^N \sum_{j=1}^N \left( \sum_{t=1}^{T-1}\log a_{ij}P(O,i_t=i,i_{t+1}=j|\overline \lambda) \right)
\]
约束条件为
\[\sum_{j=1}^Na{ij}=1
\]
引入拉格朗日乘子
\[\sum_{i=1}^N \sum_{j=1}^N \left( \sum_{t=1}^{T-1}\log a_{ij}P(O,i_t=i,i_{t+1}=j|\overline \lambda) \right)+ \gamma \left(\sum_{j=1}^Na{ij}-1\right)
\]
对\(a_{ij}\)求偏导，并令结果为0
\[\sum_{t=1}^{T-1}\frac{1}{a_{ij}}P(O,i_t=i,i_{t+1}=j|\overline \lambda)+\gamma = 0
\]
\[\sum_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\overline \lambda)+\gamma a_{ij} = 0
\]
对\(j\)求和
\[\sum_{t=1}^{T-1}P(O,i_t=i|\overline \lambda)+\gamma = 0
\]
\[\gamma = -\sum_{t=1}^{T-1}P(O,i_t=i|\overline \lambda)
\]
最后得到
\[a_{ij} = \frac{\sum\limits_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\overline \lambda)}{\sum\limits_{t=1}^{T-1}P(O,i_t=i|\overline \lambda)}
\]
- 重复2、3步直到收敛
6. 预测问题
使用维比特算法解决隐马尔可夫预测问题。基本思想：如果\(i_t^*\) 是最佳状态序列 \(I^*\)经过的状态，那么从\(i_t^*\)到\(i_T^*\)的路径必须是最优的，否则可以找出更优的\(t\)到\(T\)的路径和\(i_1\)到\(i_t^*\)拼接得到更优的路径
已知模型\(\lambda\)和观测序列\(O=\{o_1,o_2,...,o_T\}\)
\(\delta_t(i)\)表示到\(t\)时刻时（不考虑后续），状态为\(i\)的最大概率
\(\Psi_t(i)\)表示到\(t\)时刻时（不考虑后续）、使状态为\(i\)概率最大时，上一时刻的状态
\[\delta_1(i) = \pi_ib_i(o_1), \quad i = 1,2,...,N
\]
\[\delta_{t+1}(i) = \max_{1\leq j \leq N}\delta_{t}(j)a_{ji}b_i(o_t)
\]
\[\Psi_t(i) = \mathop{\arg\max}_{1\leq j \leq N} \delta_{t}(j)a_{ji}
\]
按照上述迭代至\(T\)时刻，回溯求解状态路径
\[i_T^* = \mathop{\arg\max}_{1\leq i \leq N} \Psi_T(i)
\]
\[i_{t-1}^* = \mathop{\arg\max}_{1\leq i \leq N} \Psi_t(i_t^*)
\]
最后得到观测序列和模型已知情况下概率最大的状态序列\(\{i_1^*,i_2^*,...i_T^* \}\)
import numpy as np
class HMM:
def __init__(self, state_transfer, observe, init_state, state_name, observe_name):
self.observe_name = {name: index for index,
name in enumerate(observe_name)}
self.state_name = state_name
self.A = np.array(state_transfer)
self.B = np.array(observe)
self.init_state = np.array(init_state)
self.state_nums = self.A.shape[0]
self.observe_nums = self.B.shape[1]
self.backward = None
self.forward = None
def observe_name2index(self, observe_name_list):
return [self.observe_name[name] for name in observe_name_list]
def forward_prob(self, observe_list):
ob_len = len(observe_list)
forward = np.zeros((ob_len, self.state_nums))
for i in range(self.state_nums):
forward[0][i] = self.init_state[i]*B[i][observe_list[0]]
for t in range(1, ob_len):
for i in range(self.state_nums):
current_prob = 0
for j in range(self.state_nums):
current_prob += forward[t-1][j]*self.A[j][i]
forward[t][i] = current_prob*B[i][observe_list[t]]
self.forward = forward
return forward[ob_len-1].sum()
def backward_prob(self, observe_list):
ob_len = len(observe_list)
backward = np.zeros((ob_len, self.state_nums))
for i in range(self.state_nums):
backward[ob_len-1][i] = 1
for t in range(ob_len-2, -1, -1):
for i in range(self.state_nums):
current_prob = 0
for j in range(self.state_nums):
current_prob += backward[t+1][j] * \
self.A[i][j]*self.B[j][observe_list[t+1]]
backward[t][i] = current_prob
self.backward = backward
prob = sum([self.init_state[i]*backward[0][i]*self.B[i]
[observe_list[0]] for i in range(self.state_nums)])
return prob
def observe_list_prob(self, observ_list, t=0):
if self.forward is None:
self.forward_prob(observ_list)
if self.backward is None:
self.backward_prob(observ_list)
res = 0
for i in range(self.state_nums):
for j in range(self.state_nums):
res += self.forward[t][i]*self.A[i][j] * \
self.B[j][observ_list[t+1]]*self.backward[t+1][j]
return res
def viterbi(self, observ_list):
ob_len = len(observ_list)
def func(t):
if t == 0:
prob_res = np.array([self.init_state[i]*self.B[i][observ_list[t]]
for i in range(self.state_nums)])
path_res = [[i] for i in range(self.state_nums)]
return prob_res, path_res
else:
prob_res, path_res = func(t-1)
# print("t={}:{}".format(t-1, prob_res))
next_res = []
for i in range(self.state_nums):
temp = np.array([prob_res[j]*self.A[j][i]*self.B[i][observ_list[t]]
for j in range(self.state_nums)])
index = temp.argmax()
path_res[i].append(index)
next_res.append(temp[index])
return np.array(next_res), path_res
prob_res, path_res = func(ob_len-1)
# print("t={}:{}".format(ob_len-1, prob_res))
return prob_res.argmax(), path_res[prob_res.argmax()]
if __name__ == "__main__":
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
init = [0.2, 0.4, 0.4]
Q = [1, 2, 3]
V = ['红', '白']
model = HMM(A, B, init, Q, V)
ob_list = model.observe_name2index(['红', '白', '红'])
# forward_prob = model.forward_prob(ob_list)
# backward_porb = model.backward_prob(ob_list)
# print(forward_prob, backward_porb)
# for i in range(len(ob_list)-1):
# print(model.observe_list_prob(ob_list, t=i))
prob, path = model.viterbi(ob_list)
print(prob, path)
浙公网安备 33010602011771号