机器学习—隐马尔科夫模型(7)维特比算法

776 阅读3分钟

持续创作,加速成长!这是我参与「掘金日新计划 · 6 月更文挑战」的第11天,点击查看活动详情

在给定观测序列的条件下,找到可以使观测序列出现概率最大的状态序列,也就是在观测序列和参数是已知条件下,找到一个让后验概率 P(IO,λ)P(I|O,\lambda) 最大的状态序列 II

I^=arg maxIP(IO,λ)O={o1,o2,,oT}\hat{I} = \argmax_I P(I|O,\lambda)\, O =\{o_1,o_2,\cdots,o_T\}

解码问题和预测问题

在给定观测序列的条件下,找到可以使观测序列出现概率最大的状态序列

P(i1,i2,,iTo1,o2,,oT)=P(i1,i2,,iTo1,o2,,oT)P(o1,o2,,oT)P(i_1,i_2,\cdots,i_T|o_1,o_2,\cdots,o_T) = \frac{P(i_1,i_2,\cdots,i_T|o_1,o_2,\cdots,o_T)}{P(o_1,o_2,\cdots,o_T)}

动态规划

我们将将求概率最大的问题转换为距离最短的问题,所以这里用动态规划来解决这个问题。这里要用算法就是一个比较基本动态规划算法,维特比(viterbi)算法。

δt(i)=maxi1,i2,,it1P(o1,o2,,ot,i1,i2,,it1,it=qi)\delta_t(i) = \max_{i_1,i_2,\cdots,i_{t-1}} P(o_1,o_2,\cdots,o_t,i_1,i_2,\cdots,i_{t-1},i_t=q_i)

hmm_008.png

在 t 时刻到达 qiq_i 状态的最大大概率的路径,δt\delta_t 就是表示最大值,上面公式因为 λ\lambda 是常数,所以就省略了。表示在 t 时刻前并不确定走那一条路径,只要在 t 时刻来到 qiq_i 状态即可。从而这些排列组合可能路径中选择一条概率最大一条路径来将 δt\delta_t 计算出来。

其实就是定义了局部状态,也就是在 t 时刻,当确定在 t 时刻状态为 qiq_i 所有转移状态路径中概率最大路径,也就是局部状态,要找到路径所有可能性概率最大路径。也就是我们在 t 时刻知道其状态为 qiq_i 我们是要找一条路径,这条路径到 qiq_i 的概率最大。

φt(i)=arg max1jN[δt(j)aji]\varphi_t(i) = \argmax_{1 \le j \le N} [\delta_t(j)a_{ji}]

也就是上一个节点隐状态是什么,这里 φt1(i)\varphi_{t-1}(i) 是当我们已经找到了一条路径到当前时刻 t 状态为 qiq_i 的概率最大路径,需要找到在这种情况下前一个时刻,也就是 t-1 时刻的状态是哪一个状态,这样好处是我们可以顺藤摸瓜找到 t 时刻概率最大所有对应一些列隐藏状态。

接下来根据 δt\delta_t 定义来列出 θt+1(j)\theta_{t+1}(j) 当达到 t+1 时刻状态为 qjq_j 让概率最大的路径,接下来工作就是找 δt\delta_tδt+1\delta_{t+1} 之间的关系,也就是找到递推式

δt+1(j)=maxi1,,itP(o1,o2,,ot,ot+1,i1,i2,,it1,it,it+1=qj)\delta_{t+1}(j) = \max_{i_1,\cdots,i_{t}} P(o_1,o_2,\cdots,o_t,o_{t+1},i_1,i_2,\cdots,i_{t-1},i_t,i_{t+1}=q_j)
maxδt(j)P(it+1=qjit=qi)P(ot+1it+1=qj)\max \delta_t(j) P(i_{t+1}=q_j|i_{t}=q_i)P(o_{t+1}|i_{t+1}=q_j)\\
aij=P(it+1=qjit=qi)bj(ot+1)=P(ot+1it+1=qj)a_{ij}= P(i_{t+1}=q_j|i_{t}=q_i)\\ b_j(o_{t+1}) = P(o_{t+1}|i_{t+1}=q_j)
max1iNδt(i)ajibj(ot+1)\max_{1 \le i \le N} \delta_t(i)a_{ji}b_j(o_{t+1})

也是也很好理解,也就是 j 是从 1 到 j 找到一个 δt\delta_t 最大值然后乘以从 i 到 j 转移概率在乘以状态 qjq_jot+1o_{t+1} 的生成概率作为 δt+1(j)\delta_{t+1}(j) 的值

δt(i)=P(it=qiO,λ)\delta_t(i) = P(i_t=q_i|O,\lambda)

初始化

δ1(i)=πibi(o1)i=1,2,,Nφ1(i)=0i=1,2,,N\delta_1(i) = \pi_ib_i(o_1)\,i = 1,2,\cdots,N\\ \varphi_1(i) = 0\,i = 1,2,\cdots,N
δt(i)=maxijN[δt1(j)aji]bi(ot)φt(i)=arg maxijN[δt1(j)aji]\delta_t(i) = \max_{i\le j\le N}\left[ \delta_{t-1}(j)a_{ji} \right] b_i(o_t)\\ \varphi_t(i) = \argmax_{i\le j\le N}\left[ \delta_{t-1}(j)a_{ji} \right]\\
P=maxijNδT(i)P^{*} = \max_{i\le j\le N} \delta_T(i)

这里 pp^{*} 表示最可能出现隐藏状态序列出现概率,也就是让概率最大的隐含状态

iT=arg maxijN[δT(i)]i^{*}_T = \argmax_{i\le j \le N} [\delta_T(i)]

在 T 时刻可以得到最可能隐藏状态,it=φt+1(it+1)i^*_t = \varphi_{t+1}(i^*_{t+1}) 根据 t+1t+1 时刻发生概率最大隐状态 it+1i^*_{t+1} 就可以推出前一个时刻 tt 时刻概率最大隐状态iti^*_t

实例

这里例子中,天气是状态,天气影响到 Bob 的心情,而这里观测变量是 Bob 的心情,然后在下图中列出状态转移概率和观测生成概率

hmm_007.png

为了便于计算表示我们将随机事件用具体数值表示,也就是随机变量的取值

hmm_006.png

观测序列如下,这就是我们观测 Bob 连续几天心情,由此我们来推测最近几天其所居住的城市天气状况

hmm_003.png

这里我们引入 numpy 来进行简单计算

import numpy as np

pi_0 = np.array([0.5,0.3,0.2])

A = np.array([
    [0.6,0.2,0.2],
    [0.5,0.3,0.2],
    [0.6,0.2,0.2]
])

B = np.array([
    [0.5,0.3,0.2],
    [0.3,0.3,0.4]
])

# B_0 B_1 B_0

初始化

delta_1 = pi_0*B[0]
print(delta_1) #[0.25 0.09 0.04]

viterbi_002.png

δ1(1)=π1b1(o1)=0.5×0.5=0.25δ1(1)=π1b1(o1)=0.3×0.3=0.09δ1(1)=π1b1(o1)=0.2×0.2=0.04\delta_1(1) = \pi_1b_1(o_1) = 0.5 \times 0.5 = 0.25\\ \delta_1(1) = \pi_1b_1(o_1) = 0.3 \times 0.3 = 0.09\\ \delta_1(1) = \pi_1b_1(o_1) = 0.2 \times 0.2 = 0.04\\
φ1(1)=φ1(2)=φ3(3)=0\varphi_1(1) = \varphi_1(2) = \varphi_3(3) = 0

迭代

δ2(1)=maxij3[δ1(j)aj1]b1(o2)=maxij3[0.25×0.6,0.09×0.5,0.04×0.6]×0.3=0.045δ2(2)=maxij3[δ1(j)aj1]b2(o2)=maxij3[0.25×0.2,0.09×0.3,0.04×0.2]×0.3=0.015δ2(3)=maxij3[δ1(j)aj1]b1(o2)=maxij3[0.25×0.2,0.09×0.2,0.04×0.2]×0.4=0.02\delta_2(1) = \max_{i \le j \le 3}[\delta_1(j)a_{j1}]b_1(o_2) = \max_{i \le j \le 3}[ 0.25 \times 0.6, 0.09 \times 0.5, 0.04 \times 0.6 ] \times 0.3 = 0.045\\ \delta_2(2) = \max_{i \le j \le 3}[\delta_1(j)a_{j1}]b_2(o_2) = \max_{i \le j \le 3}[ 0.25 \times 0.2, 0.09 \times 0.3, 0.04 \times 0.2 ] \times 0.3 = 0.015\\ \delta_2(3) = \max_{i \le j \le 3}[\delta_1(j)a_{j1}]b_1(o_2) = \max_{i \le j \le 3}[ 0.25 \times 0.2, 0.09 \times 0.2, 0.04 \times 0.2 ] \times 0.4 = 0.02

viterbi_005.png

delta_2 = [ np.max(delta_1 * np.transpose(A)[idx]) for idx in range(3) ]*B[1]
print(delta_2) #[0.045 0.015 0.02 ]
φ2(1)=1φ2(2)=1φ2(3)=1\varphi_2(1) = 1\\ \varphi_2(2) = 1\\ \varphi_2(3) = 1\\
print(np.argmax(delta_2))
delta_3 = [ np.max(delta_2 * np.transpose(A)[idx]) for idx in range(3) ]*B[0]
print(delta_3)#[0.0135 0.0027 0.0018]

print(np.argmax(delta_3))#0