隐马尔科夫模型三大问题

关于隐马尔科夫模型以及三大问题的相关概念，在之前一篇博客已经介绍过，这里就不再介绍了。这篇博客的主要内容是通过一个例子介绍解决隐马尔科夫模型三大问题的算法。

下面给出例子。

其中三种可观测的值是(walk,shop,clean)，隐含的状态是(rainy,sunny)。你可以简单理解为你和你的女朋友异地恋，你的女朋友在她的城市每天所做的事情和你报备，但是你不知道她的城市天气如何。

对于这个例子的三大问题描述是：
（1）已知整个模型，观测到连续三天做的事情是：walk,shop,clean，根据模型求产生这些行为的概率是多少？
（2）同样已知模型，同样连续三天所做的事情是：walk,shop,clean，求这三天的天气如何？
（3）已知这三天所做的事情是：walk,shop,clean，其他什么全都不知道，求这个模型？

遍历算法

即直接计算法，该算法用于解决第一种问题。参考李航的《统计学习方法》一书，给出此方法的数学表述。

给定模型$\lambda(A,B,\pi)$λ(A,B,π)，以及观测序列$O=(o_1,o_2,\cdots,o_T)$O=(o1​,o2​,⋯,oT​)，计算观测序列$O$O出现的概率$P(O|\lambda)$P(O∣λ)。通过列举所有可能的长度为$T$T的状态序列$I=(i_1,i_2,\cdots,i_T)$I=(i1​,i2​,⋯,iT​)，求状态序列$I$I和观测序列$O$O的联合概率$P(O,I|\lambda)$P(O,I∣λ)，然后对此求和，得到$P(O|\lambda)$P(O∣λ)。

状态序列$I=(i_1,i_2,\cdots,i_T)$I=(i1​,i2​,⋯,iT​)的概率是$P(I|\lambda)=\pi_{i_1}a_{i_1i_2}a_{i_2i_3}\cdots a_{i_{T-1}i_T}$P(I∣λ)=πi1​​ai1​i2​​ai2​i3​​⋯aiT−1​iT​​对固定的状态序列$I=(i_1,i_2,\cdots,i_T)$I=(i1​,i2​,⋯,iT​)，观测序列$O=(o_1,o_2,\cdots,o_T)$O=(o1​,o2​,⋯,oT​)的概率是$P(O|I,\lambda)=b_{i_1}(o_1)b_{i_2}(o_2)\cdots b_{i_T}(o_T)$P(O∣I,λ)=bi1​​(o1​)bi2​​(o2​)⋯biT​​(oT​)$O$O和$I$I同时出现的联合概率$P(O,I|\lambda)=P(O|I,\lambda)P(I|\lambda)$P(O,I∣λ)=P(O∣I,λ)P(I∣λ)然后求和得到$P(O|\lambda)=\sum_I{P(O,I|\lambda)}=\sum_I{P(O|I,\lambda)P(I|\lambda)}$P(O∣λ)=I∑​P(O,I∣λ)=I∑​P(O∣I,λ)P(I∣λ)

在本例中，观测序列$O$O就是walk,shop,clean（之后用w,s,c表示），而状态序列$I$I则是rainy,sunny（之后用R，S表示）的排列组合，比如说这三天可以是(R，R，R)，或者(R，S，R)。这样的排列组合总共有八种，我们只需要计算在这八种情况下(w,s,c)的概率分别是多少，然后求和就行了。

比方说，对于(R，R，R)的情况下连续三天做了(w,s,c)的概率是$P(w,s,c|R,R,R)=P(R1)*P(w|R1)*P(R2|R1)*P(s|R2)*P(R3|R1,R2)*P(c|R3)$P(w,s,c∣R,R,R)=P(R1)∗P(w∣R1)∗P(R2∣R1)∗P(s∣R2)∗P(R3∣R1,R2)∗P(c∣R3)带入数值$P(w,s,c|R,R,R)=0.6*0.1*0.7*0.4*0.7*0.5=0.0058$P(w,s,c∣R,R,R)=0.6∗0.1∗0.7∗0.4∗0.7∗0.5=0.0058这就是连续三天下雨并且做了(w,s,c)的概率。同样道理，求出其他几种情况下(w,s,c)的概率，然后求和（这里就不算了），最终得到的就是我们要求的$P(O|\lambda)$P(O∣λ)。

对于本例来说，状态序列和观测序列并不是很多，或许可以用此方法计算，但是如果状态序列和观测序列很多的情况下，该算法的计算量是很庞大的，是$O(TN^T)$O(TNT)阶的（这里$N^T$NT是指所有可能的状态序列，在本例中就是$2^3=8$23=8种可能），所以这种方法是不可行的。下面介绍计算$P(O|\lambda)$P(O∣λ)的有效算法：前向-后向算法。

前向算法

首先根据《统计学习方法》给出前向概率的定义：给定隐马尔科夫模型$\lambda$λ，定义到时刻$t$t部分观测序列为$o_1,o_2,\cdots,o_t$o1​,o2​,⋯,ot​且状态为$q_i$qi​的概率为前向概率，记为$\alpha_t(i)=P(o_1,o_2,\cdots,o_t,i_t=q_i|\lambda)$αt​(i)=P(o1​,o2​,⋯,ot​,it​=qi​∣λ)可以递推地求出前向概率$\alpha_t(i)$αt​(i)及观测序列概率$P(O|\lambda)$P(O∣λ)。

该算法可以分为以下几步：
（1）初值$\alpha_1(i)=\pi_ib_i(o_1)，i=1,2,\cdots,N$α1​(i)=πi​bi​(o1​)，i=1,2,⋯,N
（2）递推，对于$t=1,2,\cdots,T-1$t=1,2,⋯,T−1，$\alpha_{t+1}(i)=[\sum_{j=1}^N{\alpha_{t}(j)a_{ji}}]b_i(o_{t+1})，i=1,2,\cdots,N$αt+1​(i)=[j=1∑N​αt​(j)aji​]bi​(ot+1​)，i=1,2,⋯,N
（3）终止$P(O|\lambda)=\sum_{i=1}^N{\alpha_T(i)}$P(O∣λ)=i=1∑N​αT​(i)
利用前向概率计算$P(O|\lambda)$P(O∣λ)的计算量是$O(N^2T)$O(N2T)阶的，显然比直接计算的$O(TN^T)$O(TNT)要小。

接下来结合例子来做一下说明。
我们把连续三天划分为三个时刻，第一天、第二天和第三天。在第一天，我们观测到的行为是walk，它有两种可能，第一天是晴天的情况下散步和第一天是雨天的情况下散步。所以对于$t=1$t=1，计算如下：$P(w,R1)=P(R1)*P(w|R1)=0.6*0.1=0.06\\ P(w,S1)=P(S1)*P(w|S1)=0.4*0.6=0.24$P(w,R1)=P(R1)∗P(w∣R1)=0.6∗0.1=0.06P(w,S1)=P(S1)∗P(w∣S1)=0.4∗0.6=0.24
在第二天我们观测到的行为是shop，它同样有两种可能，第二天是晴天的情况下购物和第二天是雨天的情况下购物。但是第二天的行为是在第一天的基础上来的，也就是说$P(w,s,R2)=[P(w,R1)*P(R2|R1)+P(w,S1)*P(R2|S1)]*P(s|R2)\\ P(w,s,S2)=[P(w,R1)*P(S2|R1)+P(w,S1)*P(S2|S1)]*P(s|S2)\\$P(w,s,R2)=[P(w,R1)∗P(R2∣R1)+P(w,S1)∗P(R2∣S1)]∗P(s∣R2)P(w,s,S2)=[P(w,R1)∗P(S2∣R1)+P(w,S1)∗P(S2∣S1)]∗P(s∣S2)
在第三天我们观测到的行为是clean，同第二天一样，所以它的计算方式是：$P(w,s,c,R3)=[P(w,s,R2)*P(R3|R2)+P(w,s,S2)*P(R3|S2)]*P(c|R3)\\ P(w,s,c,S3)=[P(w,s,R2)*P(S3|R2)+P(w,s,S2)*P(S3|S2)]*P(c|S3)$P(w,s,c,R3)=[P(w,s,R2)∗P(R3∣R2)+P(w,s,S2)∗P(R3∣S2)]∗P(c∣R3)P(w,s,c,S3)=[P(w,s,R2)∗P(S3∣R2)+P(w,s,S2)∗P(S3∣S2)]∗P(c∣S3)
最终$P(O|\lambda)$P(O∣λ)是等于$P(w,s,c,R3)+P(w,s,c,S3)$P(w,s,c,R3)+P(w,s,c,S3)。

后向算法

依旧线给出后向概率的定义：给定隐马尔科夫模型$\lambda$λ，定义在时刻$t$t状态为$q_i$qi​的条件下，从$t+1$t+1到$T$T的部分观测序列为$o_{t+1},o_{t+2},\cdots,o_T$ot+1​,ot+2​,⋯,oT​的概率为后向概率，记为$\beta_t(i)=P(o_{t+1},o_{t+2},\cdots,o_T,i_t=q_i|\lambda)$βt​(i)=P(ot+1​,ot+2​,⋯,oT​,it​=qi​∣λ)同样可以用递推的方法求出后向概率$\beta_t(i)$βt​(i)及观测序列概率$P(O|\lambda)$P(O∣λ)。

该算法可以分为以下几步：
（1）初始化后向概率，对最终时刻的所有状态$q_i$qi​规定$\beta_T(i)=1$βT​(i)=1。
（2）对$t=T-1,T-2,\cdots,1$t=T−1,T−2,⋯,1，有$\beta_t(i)=\sum_{j=1}^N{a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}，i=1,2,\cdots,N$βt​(i)=j=1∑N​aij​bj​(ot+1​)βt+1​(j)，i=1,2,⋯,N
（3）求$P(O|\lambda)$P(O∣λ)思路与（2）一致，只不过初始概率代替了转移概率，即$P(O|\lambda)=\sum_{i=1}^N{\pi_ib_i(o_1)\beta_1(i)}$P(O∣λ)=i=1∑N​πi​bi​(o1​)β1​(i)

还是这个例子。设$\beta_3(R)=\beta_3(S)=1$β3​(R)=β3​(S)=1，开始计算$\beta_2(R)$β2​(R)和$\beta_2(S)$β2​(S)。$\beta_2(R)=a_{R->R}b_{R}(c)\beta_3(R)+a_{R->S}b_{S}(c)\beta_3(S)\\ =0.7*0.5*1+0.3*0.1*1=0.38$β2​(R)=aR−>R​bR​(c)β3​(R)+aR−>S​bS​(c)β3​(S)=0.7∗0.5∗1+0.3∗0.1∗1=0.38这个式子的意思是：假设第二天是雨天，计算第三天是雨天并且clean和第三天是晴天并且clean的情况。同理，$\beta_2(S)=a_{S->R}b_{R}(c)\beta_3(R)+a_{S->S}b_{S}(c)\beta_3(S)$β2​(S)=aS−>R​bR​(c)β3​(R)+aS−>S​bS​(c)β3​(S)
同样道理，计算$\beta_1(R)$β1​(R)和$\beta_1(S)$β1​(S)：$\beta_1(R)=a_{R->R}b_{R}(s)\beta_2(R)+a_{R->S}b_{S}(s)\beta_2(S)\\ \beta_1(S)=a_{S->R}b_{R}(s)\beta_2(R)+a_{S->S}b_{S}(s)\beta_2(S)$β1​(R)=aR−>R​bR​(s)β2​(R)+aR−>S​bS​(s)β2​(S)β1​(S)=aS−>R​bR​(s)β2​(R)+aS−>S​bS​(s)β2​(S)
最后$P(O|\lambda)=\pi_{R}b_{R}(w)\beta_1(R)+\pi_{S}b_{S}(w)\beta_1(S)$P(O∣λ)=πR​bR​(w)β1​(R)+πS​bS​(w)β1​(S)

利用前向概率和后向概率可以将观测序列概率$P(O|\lambda)$P(O∣λ)统一表示成$P(O|\lambda)=\sum_{i=1}^N{\sum_{j=1}^N{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}}，t=1,2,\cdots,T-1$P(O∣λ)=i=1∑N​j=1∑N​αt​(i)aij​bj​(ot+1​)βt+1​(j)，t=1,2,⋯,T−1此式当$t=1$t=1和$t=T-1$t=T−1时分别是两个算法第三步的两个式子。

Viterbi(维特比)算法

该算法是用来解决第二种问题。维特比算法实际上是用动态规划解隐马尔科夫模型预测问题，即用动态规划求概率最大路径（最优路径）。这时一条路径对应着一个状态序列。

根据动态规划的原理，最优路径具有这样的特性：如果最优路径在时刻$t$t通过节点$i_t^*$it∗​，那么这一路径从节点$i_t^*$it∗​到终点$i_T^*$iT∗​的部分路径，对于从$i_t^*$it∗​到$i_T^*$iT∗​的所有可能的部分路径来说，必须是最优的。因为如果不是这样，那么从$i_t^*$it∗​到$i_T^*$iT∗​就有一条更优的部分路径存在，如果把它和$i_1^*$i1∗​到终点$i_t^*$it∗​的部分路径连接起来，就会形成一条比原来路径更优的路径，这是矛盾的。

根据这一原理，我们只需要从时刻$t=1$t=1开始，递推地计算在时刻$t$t状态为$i$i的各条部分路径的概率最大值，直至得到时刻$t=T$t=T状态为$i$i的各条部分路径的最大概率，时刻$t=T$t=T的最大概率即为最优路径的概率$P^*$P∗，最优路径终点$i_T^*$iT∗​也同时得到。之后，为了找出最优路径的各个节点，从终点$i_T^*$iT∗​开始，向前逐步求得$i_{T-1}^*,\cdots,i_1^*$iT−1∗​,⋯,i1∗​，得到最优路径$I^*=(i_1^*,i_2^*,\cdots,i_T^*)$I∗=(i1∗​,i2∗​,⋯,iT∗​)。这就是维特比算法。

我们引入两个变量$\Delta$Δ和$\phi$ϕ。定义在时刻$t$t状态为$i$i的所有单个路径$(i_1,i_2,\cdots,i_t)$(i1​,i2​,⋯,it​)中的概率最大值为$\Delta_t(i)=\mathop{\text{max}}\limits_{i_1,i_2,\cdots,i_{t-1}}P(i_t=i,i_{t-1},\cdots,i_1,o_t,\cdots,o_1|\lambda)，i=1,2,\cdots,N$Δt​(i)=i1​,i2​,⋯,it−1​max​P(it​=i,it−1​,⋯,i1​,ot​,⋯,o1​∣λ)，i=1,2,⋯,N由定义可得变量$\sigma$σ的递推公式$\Delta_{t+1}(i)=\mathop{\text{max}}\limits_{i_1,i_2,\cdots,i_{t}}P(i_{t+1}=i,i_{t},\cdots,i_1,o_{t+1},\cdots,o_1|\lambda)\\ =\mathop{\text{max}}\limits_{1\leq j \leq N}[\Delta_{t}(j)a_{ji}]b_i(o_{t+1})，i=1,2,\cdots,N；t=1,2,\cdots,T-1$Δt+1​(i)=i1​,i2​,⋯,it​max​P(it+1​=i,it​,⋯,i1​,ot+1​,⋯,o1​∣λ)=1≤j≤Nmax​[Δt​(j)aji​]bi​(ot+1​)，i=1,2,⋯,N；t=1,2,⋯,T−1定义在时刻$t$t状态为$i$i的所有单个路径$(i_1,i_2,\cdots,i_t)$(i1​,i2​,⋯,it​)中概率最大的路径的第$i-1$i−1个节点为$\phi_t(i)=\mathop{\text{argmax}}\limits_{1\leq j \leq N}[\Delta_{t-1}(j)a_{ji}]，i=1,2,\cdots,N$ϕt​(i)=1≤j≤Nargmax​[Δt−1​(j)aji​]，i=1,2,⋯,N

维特比算法的步骤如下：
（1）初始化$\Delta_1(i)=\pi_ib_i(o_1)，i=1,2,\cdots,N\\ \phi_1(i)=0，i=1,2,\cdots,N$Δ1​(i)=πi​bi​(o1​)，i=1,2,⋯,Nϕ1​(i)=0，i=1,2,⋯,N
（2）递推。对于$t=2,3,\cdots,T$t=2,3,⋯,T，$\Delta_{t}(i)=\mathop{\text{max}}\limits_{1\leq j \leq N}[\Delta_{t-1}(j)a_{ji}]b_i(o_{t})，i=1,2,\cdots,N\\ \phi_t(i)=\mathop{\text{argmax}}\limits_{1\leq j \leq N}[\Delta_{t-1}(j)a_{ji}]，i=1,2,\cdots,N$Δt​(i)=1≤j≤Nmax​[Δt−1​(j)aji​]bi​(ot​)，i=1,2,⋯,Nϕt​(i)=1≤j≤Nargmax​[Δt−1​(j)aji​]，i=1,2,⋯,N
（3）终止$P^*=\mathop{\text{max}}\limits_{1\leq i \leq N}\Delta_{T}(i)\\ i_T^*=\mathop{\text{argmax}}\limits_{1\leq i \leq N}[\Delta_{T}(i)]$P∗=1≤i≤Nmax​ΔT​(i)iT∗​=1≤i≤Nargmax​[ΔT​(i)]
（4）最优路径回溯。对于$t=T-1,T-2,\cdots,1$t=T−1,T−2,⋯,1，$i_t^*=\phi_{t+1}(i_{t+1}^*)$it∗​=ϕt+1​(it+1∗​)于是，求得最优路径$I^*=(i_1^*,i_2^*,\cdots,i_T^*)$I∗=(i1∗​,i2∗​,⋯,iT∗​)。

对于本例，我们设$\phi_1(R)=\phi_1(S)=0$ϕ1​(R)=ϕ1​(S)=0，这是因为对于第一天的晴天或雨天状态来说它没有前驱。然后计算$\Delta_1(R)$Δ1​(R)和$\Delta_1(S)$Δ1​(S)，这两个变量表示的是到目前为止最高可能性是多少，就是说第一天是雨天的可能性是多少和第一天是晴天的可能性是多少。$\Delta_1(R)=\pi_R*b_R(w)=0.6*0.1=0.06\\ \Delta_1(S)=\pi_S*b_S(w)=0.4*0.6=0.24$Δ1​(R)=πR​∗bR​(w)=0.6∗0.1=0.06Δ1​(S)=πS​∗bS​(w)=0.4∗0.6=0.24
接下来我们计算$\Delta_2(R)$Δ2​(R)和$\Delta_2(S)$Δ2​(S)，$\Delta_2(R)=max(\Delta_1(R)*a_{R->R},\Delta_1(S)*a_{S->R})*b_R(s)\\ =max(0.06*0.7,0.24*0.4)*0.4=0.0384$Δ2​(R)=max(Δ1​(R)∗aR−>R​,Δ1​(S)∗aS−>R​)∗bR​(s)=max(0.06∗0.7,0.24∗0.4)∗0.4=0.0384因为$\Delta_1(R)*a_{R->R}$Δ1​(R)∗aR−>R​比$\Delta_1(S)*a_{S->R}$Δ1​(S)∗aS−>R​要小，所以$\phi_2(R)=S$ϕ2​(R)=S，就是说在第二天是雨天的情况下它的前序最优解是晴天。同样的方法求出$\Delta_2(S)$Δ2​(S)、$\phi_2(S)$ϕ2​(S)、$\Delta_3(R)$Δ3​(R)、$\phi_3(R)$ϕ3​(R)、$\Delta_3(S)$Δ3​(S)、$\phi_3(S)$ϕ3​(S)。

最后我们比较$\Delta_3(R)$Δ3​(R)和$\Delta_3(S)$Δ3​(S)。假设$\Delta_3(R)>\Delta_3(S)$Δ3​(R)>Δ3​(S)，那么可以认为连续三天做了(w,s,c)的最优解在第三天是雨天。往上回溯，假设$\phi_3(R)=S$ϕ3​(R)=S（这里是因为上面没做计算），那么最优解在第二天是晴天，继续往上，假设$\phi_2(S)=R$ϕ2​(S)=R（同样没计算），那么最优解在第一天是雨天，而$\phi_1(R)=0$ϕ1​(R)=0，所以回溯结束。最终按照倒序排列得到$(R,S,R)$(R,S,R)，即第一天是雨天第二天是晴天第三天是雨天，这就是所有的最优解。

Baum-Welch算法

这个算法是用来解决第三种问题。假定给定的训练数据只包含$S$S个长度为$T$T的观测序列$O_1,O_2,\cdots,O_S$O1​,O2​,⋯,OS​而没有对应的状态序列，目标是学习隐马尔科夫模型$\lambda=(A,B,\pi)$λ=(A,B,π)的参数。我们把观测数据看做观测数据$O$O，状态序列数据看做不可观测的隐藏数据$I$I，那么隐马尔科夫模型实际上就是一个含有隐藏变量的概率模型$P(O|\lambda)=\sum_I{P(O|I,\lambda)P(I|\lambda)}$P(O∣λ)=I∑​P(O∣I,λ)P(I∣λ)它的参数学习可以由EM算法实现。无奈没有学过EM算法，所以只好继续挖坑，下次再填吧。