1. HMM模型参数求解概述
HMM模型参数求解根据已知的条件可以分为两种情况。
第一种情况较为简单，就是我们已知 D D D个长度为 T T T的观测序列和对应的隐藏状态序列，即 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , . . . ( O D , I D ) } \{(O_1, I_1), (O_2, I_2), ...(O_D, I_D)\} {(O1,I1),(O2,I2),...(OD,ID)}是已知的，此时我们可以很容易的用最大似然来求解模型参数。
假设样本从隐藏状态 q i q_i qi转移到 q j q_j qj的频率计数是 A i j A_{ij} Aij,那么状态转移矩阵求得为： A = [ a i j ] , 其 中 a i j = A i j ∑ s = 1 N A i s A = \Big[a_{ij}\Big], \;其中a_{ij} = \frac{A_{ij}}{\sum\limits_{s=1}^{N}A_{is}} A=[aij],其中aij=s=1∑NAisAij
假设样本隐藏状态为 q j q_j qj且观测状态为 v k v_k vk的频率计数是 B j k B_{jk} Bjk,那么观测状态概率矩阵为： B = [ b j ( k ) ] , 其 中 b j ( k ) = B j k ∑ s = 1 M B j s B= \Big[b_{j}(k)\Big], \;其中b_{j}(k)= \frac{B_{jk}}{\sum\limits_{s=1}^{M}B_{js}} B=[bj(k)],其中bj(k)=s=1∑MBjsBjk
假设所有样本中初始隐藏状态为 q i q_i qi的频率计数为 C ( i ) C(i) C(i),那么初始概率分布为： Π = π ( i ) = C ( i ) ∑ s = 1 N C ( s ) \Pi = \pi(i) = \frac{C(i)}{\sum\limits_{s=1}^{N}C(s)} Π=π(i)=s=1∑NC(s)C(i)
可见第一种情况下求解模型还是很简单的。但是在很多时候，我们无法得到HMM样本观察序列对应的隐藏序列，只有 D D D个长度为 T T T的观测序列，即 { ( O 1 ) , ( O 2 ) , . . . ( O D ) } \{(O_1), (O_2), ...(O_D)\} {(O1),(O2),...(OD)}是已知的，此时我们能不能求出合适的HMM模型参数呢？这就是我们的第二种情况，也是我们本文要讨论的重点。它的解法最常用的是鲍姆-韦尔奇算法，其实就是基于EM算法的求解，只不过鲍姆-韦尔奇算法出现的时代，EM算法还没有被抽象出来，所以我们本文还是说鲍姆-韦尔奇算法。
2. 鲍姆-韦尔奇算法原理
鲍姆-韦尔奇算法原理既然使用的就是EM算法的原理，那么我们需要在E步求出联合分布 P ( O , I ∣ λ ) P(O,I|\lambda) P(O,I∣λ)基于条件概率 P ( I ∣ O , λ ‾ ) P(I|O,\overline{\lambda}) P(I∣O,λ)的期望，其中 λ ‾ \overline{\lambda} λ为当前的模型参数，然后再M步最大化这个期望，得到更新的模型参数 λ \lambda λ。接着不停的进行EM迭代，直到模型参数的值收敛为止。
首先来看看E步，当前模型参数为 λ ‾ \overline{\lambda} λ, 联合分布 P ( O , I ∣ λ ) P(O,I|\lambda) P(O,I∣λ)基于条件概率 P ( I ∣ O , λ ‾ ) P(I|O,\overline{\lambda}) P(I∣O,λ)的期望表达式为： L ( λ , λ ‾ ) = ∑ I P ( I ∣ O , λ ‾ ) l o g P ( O , I ∣ λ ) L(\lambda, \overline{\lambda}) = \sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda) L(λ,λ)=I∑P(I∣O,λ)logP(O,I∣λ)
在M步，我们极大化上式，然后得到更新后的模型参数如下： λ ‾ = a r g max λ ∑ I P ( I ∣ O , λ ‾ ) l o g P ( O , I ∣ λ ) \overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda) λ=argλmaxI∑P(I∣O,λ)logP(O,I∣λ)
通过不断的E步和M步的迭代，直到 λ ‾ \overline{\lambda} λ收敛。下面我们来看看鲍姆-韦尔奇算法的推导过程。
3. 鲍姆-韦尔奇算法的推导
我们的训练数据为 { ( O 1 , I 1 ) , ( O 2 , I 2 ) , . . . ( O D , I D ) } \{(O_1, I_1), (O_2, I_2), ...(O_D, I_D)\} {(O1,I1),(O2,I2),...(OD,ID)}，其中任意一个观测序列 O d = { o 1 ( d ) , o 2 ( d ) , . . . o T ( d ) } O_d = \{o_1^{(d)}, o_2^{(d)}, ... o_T^{(d)}\} Od={o1(d),o2(d),...oT(d)},其对应的未知的隐藏状态序列表示为： I d = { i 1 ( d ) , i 2 ( d ) , . . . i T ( d ) } I_d = \{i_1^{(d)}, i_2^{(d)}, ... i_T^{(d)}\} Id={i1(d),i2(d),...iT(d)}
首先看鲍姆-韦尔奇算法的E步，我们需要先计算联合分布 P ( O , I ∣ λ ) P(O,I|\lambda) P(O,I∣λ)的表达式如下： P ( O , I ∣ λ ) = ∏ d = 1 D π i 1 ( d ) b i 1 ( d ) ( o 1 ( d ) ) a i 1 ( d ) i 2 ( d ) b i 2 ( d ) ( o 2 ( d ) ) . . . a i T − 1 ( d ) i T ( d ) b i T ( d ) ( o T ( d ) ) P(O,I|\lambda) = \prod_{d=1}^D\pi_{i_1^{(d)}}b_{i_1^{(d)}}(o_1^{(d)})a_{i_1^{(d)}i_2^{(d)}}b_{i_2^{(d)}}(o_2^{(d)})...a_{i_{T-1}^{(d)}i_T^{(d)}}b_{i_T^{(d)}}(o_T^{(d)}) P(O,I∣λ)=d=1∏Dπi1(d)bi1(d)(o1(d))ai1(d)i2(d)bi2(d)(o2(d))...aiT−1(d)iT(d)biT(d)(oT(d))
我们的E步得到的期望表达式为： L ( λ , λ ‾ ) = ∑ I P ( I ∣ O , λ ‾ ) l o g P ( O , I ∣ λ ) L(\lambda,\overline{\lambda}) = \sum\limits_{I}P(I|O,\overline{\lambda})logP(O,I|\lambda) L(λ,λ)=I∑P(I∣O,λ)logP(O,I∣λ)
在M步我们要极大化上式。由于 P ( I ∣ O , λ ‾ ) = P ( I , O ∣ λ ‾ ) / P ( O ∣ λ ‾ ) P(I|O,\overline{\lambda}) = P(I,O|\overline{\lambda})/P(O|\overline{\lambda}) P(I∣O,λ)=P(I,O∣λ)/P(O∣λ),而 P ( O ∣ λ ‾ ) P(O|\overline{\lambda}) P(O∣λ)是常数，因此我们要极大化的式子等价于： λ ‾ = a r g max λ ∑ I P ( O , I ∣ λ ‾ ) l o g P ( O , I ∣ λ ) \overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{I}P(O,I|\overline{\lambda})logP(O,I|\lambda) λ=argλmaxI∑P(O,I∣λ)logP(O,I∣λ)
我们将上面 P ( O , I ∣ λ ) P(O,I|\lambda) P(O,I∣λ)的表达式带入我们的极大化式子，得到的表达式如下： λ ‾ = a r g max λ ∑ d = 1 D ∑ I P ( O , I ∣ λ ‾ ) ( l o g π i 1 + ∑ t = 1 T − 1 l o g a i t , i t + 1 + ∑ t = 1 T l o g b i t ( o t ) ) \overline{\lambda} = arg\;\max_{\lambda}\sum\limits_{d=1}^D\sum\limits_{I}P(O,I|\overline{\lambda})(log\pi_{i_1} + \sum\limits_{t=1}^{T-1}log\;a_{i_t,i_{t+1}} + \sum\limits_{t=1}^Tlog b_{i_t}(o_t)) λ=argλmaxd=1∑DI∑P(O,I∣λ)(logπi1+t=1∑T−1logait,it+1+t=1∑Tlogbit(ot))
我们的隐藏模型参数 λ = ( A , B , Π ) \lambda =(A,B,\Pi) λ=(A,B,Π),因此下面我们只需要对上式分别对 A , B , Π A,B,\Pi A,B,Π求导即可得到我们更新的模型参数 λ ‾ \overline{\lambda} λ
首先我们看看对模型参数 Π \Pi Π的求导。由于 Π \Pi Π只在上式中括号里的第一部分出现，因此我们对于 Π \Pi Π的极大化式子为： π i ‾ = a r g max π i 1 ∑ d = 1 D ∑ I P ( O , I ∣ λ ‾ ) l o g π i 1 = a r g max π i ∑ d = 1 D ∑ i = 1 N P ( O , i 1 ( d ) = i ∣ λ ‾ ) l o g π i \overline{\pi_i} = arg\;\max_{\pi_{i_1}} \sum\limits_{d=1}^D\sum\limits_{I}P(O,I|\overline{\lambda})log\pi_{i_1} = arg\;\max_{\pi_{i}} \sum\limits_{d=1}^D\sum\limits_{i=1}^NP(O,i_1^{(d)} =i|\overline{\lambda})log\pi_{i} πi=argπi1maxd=1∑DI∑P(O,I∣λ)logπi1=argπimaxd=1∑Di=1∑NP(O,i1(d)=i∣λ)logπi
由于 π i \pi_i πi还满足 ∑ i = 1 N π i = 1 \sum\limits_{i=1}^N\pi_i =1 i=1∑Nπi=1，因此根据拉格朗日子乘法，我们得到 π i \pi_i πi要极大化的拉格朗日函数为： a r g max π i ∑ d = 1 D ∑ i = 1 N P ( O , i 1 ( d ) = i ∣ λ ‾ ) l o g π i + γ ( ∑ i = 1 N π i − 1 ) arg\;\max_{\pi_{i}}\sum\limits_{d=1}^D\sum\limits_{i=1}^NP(O,i_1^{(d)} =i|\overline{\lambda})log\pi_{i} + \gamma(\sum\limits_{i=1}^N\pi_i -1) argπimaxd=1∑Di=1∑NP(O,i1(d)=i∣λ)logπi+γ(i=1∑Nπi−1)
其中， γ \gamma γ为拉格朗日系数。上式对 π i \pi_i πi求偏导数并令结果为0， 我们得到： ∑ d = 1 D P ( O , i 1 ( d ) = i ∣ λ ‾ ) + γ π i = 0 \sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda}) + \gamma\pi_i = 0 d=1∑DP(O,i1(d)=i∣λ)+γπi=0
令 i i i分别等于从1到 N N N，从上式可以得到 N N N个式子，对这 N N N个式子求和可得： ∑ d = 1 D P ( O ∣ λ ‾ ) + γ = 0 \sum\limits_{d=1}^DP(O|\overline{\lambda}) + \gamma = 0 d=1∑DP(O∣λ)+γ=0
从上两式消去 γ \gamma γ,得到 π i \pi_i πi的表达式为： π i = ∑ d = 1 D P ( O , i 1 ( d ) = i ∣ λ ‾ ) ∑ d = 1 D P ( O ∣ λ ‾ ) = ∑ d = 1 D P ( O , i 1 ( d ) = i ∣ λ ‾ ) D P ( O ∣ λ ‾ ) = ∑ d = 1 D P ( i 1 ( d ) = i ∣ O , λ ‾ ) D = ∑ d = 1 D P ( i 1 ( d ) = i ∣ O ( d ) , λ ‾ ) D \pi_i =\frac{\sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda})}{\sum\limits_{d=1}^DP(O|\overline{\lambda})} = \frac{\sum\limits_{d=1}^DP(O,i_1^{(d)} =i|\overline{\lambda})}{DP(O|\overline{\lambda})} = \frac{\sum\limits_{d=1}^DP(i_1^{(d)} =i|O, \overline{\lambda})}{D} = \frac{\sum\limits_{d=1}^DP(i_1^{(d)} =i|O^{(d)}, \overline{\lambda})}{D} πi=d=1∑DP(O∣λ)d=1∑DP(O,i1(d)=i∣λ)=DP(O∣λ)d=1∑DP(O,i1(d)=i∣λ)=Dd=1∑DP(i1(d)=i∣O,λ)=Dd=1∑DP(i1(d)=i∣O(d),λ)
利用我们在隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义可得： P ( i 1 ( d ) = i ∣ O ( d ) , λ ‾ ) = γ 1 ( d ) ( i ) P(i_1^{(d)} =i|O^{(d)}, \overline{\lambda}) = \gamma_1^{(d)}(i) P(i1(d)=i∣O(d),λ)=γ1(d)(i)
因此最终我们在M步 π i \pi_i πi的迭代公式为： π i = ∑ d = 1 D γ 1 ( d ) ( i ) D \pi_i = \frac{\sum\limits_{d=1}^D\gamma_1^{(d)}(i)}{D} πi=Dd=1∑Dγ1(d)(i)
现在我们来看看 A A A的迭代公式求法。方法和 Π \Pi Π的类似。由于 A A A只在最大化函数式中括号里的第二部分出现，而这部分式子可以整理为： ∑ d = 1 D ∑ I ∑ t = 1 T − 1 P ( O , I ∣ λ ‾ ) l o g a i t , i t + 1 = ∑ d = 1 D ∑ i = 1 N ∑ j = 1 N ∑ t = 1 T − 1 P ( O , i t ( d ) = i , i t + 1 ( d ) = j ∣ λ ‾ ) l o g a i j \sum\limits_{d=1}^D\sum\limits_{I}\sum\limits_{t=1}^{T-1}P(O,I|\overline{\lambda})log\;a_{i_t,i_{t+1}} = \sum\limits_{d=1}^D\sum\limits_{i=1}^N\sum\limits_{j=1}^N\sum\limits_{t=1}^{T-1}P(O,i_t^{(d)} = i, i_{t+1}^{(d)} = j|\overline{\lambda})log\;a_{ij} d=1∑DI∑t=1∑T−1P(O,I∣λ)logait,it+1=d=1∑Di=1∑Nj=1∑Nt=1∑T−1P(O,it(d)=i,it+1(d)=j∣λ)logaij
由于 a i j a_{ij} aij还满足 ∑ j = 1 N a i j = 1 \sum\limits_{j=1}^Na_{ij} =1 j=1∑Naij=1。和求解 π i \pi_i πi类似，我们可以用拉格朗日子乘法并对 a i j a_{ij} aij求导，并令结果为0，可以得到 a i j a_{ij} aij的迭代表达式为： a i j = ∑ d = 1 D ∑ t = 1 T − 1 P ( O ( d ) , i t ( d ) = i , i t + 1 ( d ) = j ∣ λ ‾ ) ∑ d = 1 D ∑ t = 1 T − 1 P ( O ( d ) , i t ( d ) = i ∣ λ ‾ ) a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}P(O^{(d)}, i_t^{(d)} = i, i_{t+1}^{(d)} = j|\overline{\lambda})}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}P(O^{(d)}, i_t^{(d)} = i|\overline{\lambda})} aij=d=1∑Dt=1∑T−1P(O(d),it(d)=i∣λ)d=1∑Dt=1∑T−1P(O(d),it(d)=i,it+1(d)=j∣λ)
利用隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义和第五节 ξ t ( i , j ) \xi_t(i,j) ξt(i,j)的定义可得们在M步 a i j a_{ij} aij的迭代公式为： a i j = ∑ d = 1 D ∑ t = 1 T − 1 ξ t ( d ) ( i , j ) ∑ d = 1 D ∑ t = 1 T − 1 γ t ( d ) ( i ) a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\xi_t^{(d)}(i,j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\gamma_t^{(d)}(i)} aij=d=1∑Dt=1∑T−1γt(d)(i)d=1∑Dt=1∑T−1ξt(d)(i,j)
现在我们来看看 B B B的迭代公式求法。方法和 Π \Pi Π的类似。由于 B B B只在最大化函数式中括号里的第三部分出现，而这部分式子可以整理为： ∑ d = 1 D ∑ I ∑ t = 1 T P ( O , I ∣ λ ‾ ) l o g b i t ( o t ) = ∑ d = 1 D ∑ j = 1 N ∑ t = 1 T P ( O , i t ( d ) = j ∣ λ ‾ ) l o g b j ( o t ) \sum\limits_{d=1}^D\sum\limits_{I}\sum\limits_{t=1}^{T}P(O,I|\overline{\lambda})log\;b_{i_t}(o_t) = \sum\limits_{d=1}^D\sum\limits_{j=1}^N\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})log\;b_{j}(o_t) d=1∑DI∑t=1∑TP(O,I∣λ)logbit(ot)=d=1∑Dj=1∑Nt=1∑TP(O,it(d)=j∣λ)logbj(ot)
由于 b j ( o t ) b_{j}(o_t) bj(ot)还满足 ∑ k = 1 M b j ( o t = v k ) = 1 \sum\limits_{k=1}^Mb_{j}(o_t =v_k) =1 k=1∑Mbj(ot=vk)=1。和求解 π i \pi_i πi类似，我们可以用拉格朗日子乘法并对 b j ( k ) b_{j}(k) bj(k)求导，并令结果为0，得到 b j ( k ) b_{j}(k) bj(k)的迭代表达式为： b j ( k ) = ∑ d = 1 D ∑ t = 1 T P ( O , i t ( d ) = j ∣ λ ‾ ) I ( o t ( d ) = v k ) ∑ d = 1 D ∑ t = 1 T P ( O , i t ( d ) = j ∣ λ ‾ ) b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})I(o_t^{(d)}=v_k)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}P(O,i_t^{(d)} = j|\overline{\lambda})} bj(k)=d=1∑Dt=1∑TP(O,it(d)=j∣λ)d=1∑Dt=1∑TP(O,it(d)=j∣λ)I(ot(d)=vk)
其中 I ( o t ( d ) = v k ) I(o_t^{(d)}=v_k) I(ot(d)=vk)当且仅当 o t ( d ) = v k o_t^{(d)}=v_k ot(d)=vk时为1，否则为0. 利用隐马尔科夫模型HMM（二）前向后向算法评估观察序列概率里第二节中前向概率的定义可得 b j ( o t ) b_{j}(o_t) bj(ot)的最终表达式为： b j ( k ) = ∑ d = 1 D ∑ t = 1 , o t ( d ) = v k T γ t ( d ) ( j ) ∑ d = 1 D ∑ t = 1 T γ t ( d ) ( j ) b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1, o_t^{(d)}=v_k}^{T}\gamma_t^{(d)}(j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}\gamma_t^{(d)}(j)} bj(k)=d=1∑Dt=1∑Tγt(d)(j)d=1∑Dt=1,ot(d)=vk∑Tγt(d)(j)
有了 π i , a i j , b j ( k ) \pi_i, a_{ij},b_{j}(k) πi,aij,bj(k)的迭代公式，我们就可以迭代求解HMM模型参数了。
4. 鲍姆-韦尔奇算法流程总结
这里我们概括总结下鲍姆-韦尔奇算法的流程。
输入： D D D个观测序列样本 { ( O 1 ) , ( O 2 ) , . . . ( O D ) } \{(O_1), (O_2), ...(O_D)\} {(O1),(O2),...(OD)}
输出：HMM模型参数
1)随机初始化所有的 π i , a i j , b j ( k ) \pi_i, a_{ij},b_{j}(k) πi,aij,bj(k)
2) 对于每个样本 d = 1 , 2 , . . . D d = 1,2,...D d=1,2,...D，用前向后向算法计算 γ t ( d ) ( i ) ， ξ t ( d ) ( i , j ) , t = 1 , 2... T \gamma_t^{(d)}(i)，\xi_t^{(d)}(i,j), t =1,2...T γt(d)(i)，ξt(d)(i,j),t=1,2...T
3) 更新模型参数：
π
i
=
∑
d
=
1
D
γ
1
(
d
)
(
i
)
D
\pi_i = \frac{\sum\limits_{d=1}^D\gamma_1^{(d)}(i)}{D}
πi=Dd=1∑Dγ1(d)(i)
a
i
j
=
∑
d
=
1
D
∑
t
=
1
T
−
1
ξ
t
(
d
)
(
i
,
j
)
∑
d
=
1
D
∑
t
=
1
T
−
1
γ
t
(
d
)
(
i
)
a_{ij} = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\xi_t^{(d)}(i,j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T-1}\gamma_t^{(d)}(i)}
aij=d=1∑Dt=1∑T−1γt(d)(i)d=1∑Dt=1∑T−1ξt(d)(i,j)
b
j
(
k
)
=
∑
d
=
1
D
∑
t
=
1
,
o
t
(
d
)
=
v
k
T
γ
t
(
d
)
(
j
)
∑
d
=
1
D
∑
t
=
1
T
γ
t
(
d
)
(
j
)
b_{j}(k) = \frac{\sum\limits_{d=1}^D\sum\limits_{t=1, o_t^{(d)}=v_k}^{T}\gamma_t^{(d)}(j)}{\sum\limits_{d=1}^D\sum\limits_{t=1}^{T}\gamma_t^{(d)}(j)}
bj(k)=d=1∑Dt=1∑Tγt(d)(j)d=1∑Dt=1,ot(d)=vk∑Tγt(d)(j)
4) 如果 π i , a i j , b j ( k ) \pi_i, a_{ij},b_{j}(k) πi,aij,bj(k)的值已经收敛，则算法结束，否则回到第2步继续迭代。
以上就是鲍姆-韦尔奇算法的整个过程。
（文章转自:https://www.cnblogs.com/pinard/p/6972299.html）