HMM隐马尔科夫模型

①通俗的理解

首先举一个例子，扔骰子，有三种骰子，第一个是比较常见的6个面x = [1,2,3,4,5,6]，每一个面的概率都是1/6。第二个只有4个面，x = [1,2,3,4]，每一个面的概率是1/4。第三个有8个面，x = [1,2,3,4,5,6,7,8]，每一个面的概率是1/8。

首先先选择一个骰子，挑到每一个骰子的概率1/3，然后投掷，可以得到x = [1,2,3,4,5,6,7,8]。最后会得到一堆的序列，比如会得到O = [1,5,3,8,6,5,7,2]等等，这种序列叫可见状态序列，但在HMM里面，还存在一个隐含状态链，比如这个状态链可能是I = [D6,D8,D4,D6,D8]。从8个面出来的，6个面投出来的。所以隐马尔科夫模型里面有两串序列，可观测序列和隐含状态序列。
一般来讲，说到隐马尔科夫链其实就是隐含的状态序列了，markov链各个元素直接是存在了转化关系的，比如一开始是D6下一个状态是D4，D6，D8的概率都是1/3，这样设置平均值只是为了初始化而已，一般这种值其实是可以随意改变的，比如这里下一个状态序列是D8，那么就可以设置转化到D6的概率是0.1，D4的概率是0.1，D8的概率是0.8。这样就又是一个很新的HMM了。
同样的，尽管可见的状态之间没有转换概率，但是隐含状态和可见状态之间有一个概率叫做输出概率，也叫发射概率，就这个例子来说，D6输出1的概率就是6了，产生23456的概率都是1/6。
所以，HMM首先是有两个序列，可观测序列，隐含序列。每一个隐含序列里面的元素直接是存在着转化概率的，也就是用概率来描述转化到下一个状态的可能性；而隐含序列和观测序列之间也有一个概率，发射概率，也叫输出概率。

隐含状态关系图：

每一个箭头都会对应着一个概率，转换概率。
如果知道了转换概率和隐含概率那么求解观测序列会很简单，但如果是这样那么这个模型就没有什么可以研究的了，所以肯定是会缺失一点值的，比如知道观测序列和参数，但是不知道隐含序列是什么，或者是知道了 观测序列想求参数，又或者是知道了这个参数和隐含序列求这个观测序列出现的参数是什么。其实这些就对应了后面的三个问题。HMM的核心其实也就是这三个问题。

②三个基本问题通俗解释

1.知道骰子有几种（知道隐含序列有多少种），也就是几种骰子，上面就有三种，每种骰子是什么（转换概率），根据骰子扔出的结果（可观测状态链），想知道每次投的是哪种骰子（隐含状态序列）？
这问题其实就是解码问题了，知道观测序列求隐含序列。其实有两种解法，一种就是暴力求解：其实就是maximum likelihood全部乘起来一起算即可。另一种就厉害带了，递推的方法，前向算法，后向算法，前向-后向算法。
2.还是知道骰子有几种隐含状态（隐含状态的数量），知道骰子是什么（知道转换概率），根据骰子投出的结果我想知道投出这个结果的概率是多少？
这个问题看起来貌似没有什么太大的作用，但是实际上是用于对于模型的检测，如果我们投出来的结果都对应了很小的概率，那么模型可能就是错误的，有人换了我们的骰子。
3.知道骰子有几种（知道隐含序列的种类），不知道每种骰子是什么（不知道转换概率），实验得到多次骰子的结果（可观测序列），现在想知道每种骰子是什么（转换概率）
这很重要，后面会用到来求解分词的状态。
到这里先引出一个比较简单的问题，0号问题：
0.知道骰子的种类，骰子是什么，每次扔什么骰子，根据这个结果，求投出这个结果的概率是多少。
其实就是该知道的都知道了，直接求概率。

求解无非就是概率相乘了：
P = P(D6)*P(D6->1)*P(D6->D8)*P(D8->6)*P(D8-D8)*P(D8->3) = \frac{1}{3}*\frac{1}{6}*\frac{1}{3}*\frac{1}{8}*\frac{1}{3}*\frac{1}{8}

③三个问题的简单解答

1.看见了观测序列，求隐藏序列

这里解释的就是第一种解法了，最大似然估计，也就是暴力解法。首先我知道我有3个骰子，六面骰子，八面骰子，四面骰子，同时结果我也知道了（1 6 3 5 2 7 3 5 2 4），现在想知道我每一次投的哪一个骰子，是六面的还是四面的还是八面的呢？
最简单的方法就是穷举了，从零个一直到第最后一个一次把每一个概率算出来，链不长还行，链要是长了就算不了了，穷举的数量太大。
接下来是要讨论另一个比较牛逼的算法，viterbi algorithm。
首先只扔一次骰子：

结果是1的话那么概率最大的就是四面骰子了，1/4，其他的都是1/6和1/8。接着再扔一次： 这个时候我们就要计算三个值了，分别是D6，D8，D4的概率。要取到最大的概率，那么第一次肯定就是取到D4了，所以取到D6的最大概率：P_2(D6) = P(D4)*P(D4->1)*P(D6)*P(D6->6) = \frac{1}{3}*\frac{1}{4}*\frac{1}{3}*\frac{1}{6}
取到D8的最大概率：P_2(D8) = P(D4)*P(D4->1)*P(D8)*P(D8->6) = \frac{1}{3}*\frac{1}{4}*\frac{1}{3}*\frac{1}{8}取到D4这辈子都不可能了，四面骰子没有6。所以最大的序列就是D4，D6
继续拓展一下，投多一个： 这个时候又要计算三种情况了，分别计算为D4，D6，D8的概率是多少：
P_3(D4) = P_2(D6)*P(D6->D4)*P(D4->3) = \frac{1}{216}*\frac{1}{3}*\frac{1}{4}
P_3(D6) = P_2(D6)*P(D6->D6)*P(D6->3) = \frac{1}{216}*\frac{1}{3}*\frac{1}{6}
P_3(D8) = P_2(D6)*P(D6->D8)*P(D8->3) = \frac{1}{216}*\frac{1}{3}*\frac{1}{8}
可以看到，概率最大的就是D4了，所以三次投骰子概率最大的隐含序列就是D4，D6，D4。所以无论这个序列多长，直接递推计算可以算出来的，算到最后一位的时候，直接反着找就好了。

2.求可观测序列的概率

如果你的骰子有问题，要怎么检测出来？首先，可以先算一下正常骰子投出一段序列的概率，再算算不正常的投出来序列的概率即可。如果小于正常的，那么就有可能是错的了。
比如结果如上，我们这个时候就不是要求隐含序列了，而是求出现这个观测序列的总概率是多少。首先如果是只投一个骰子： 这时候三种骰子都有可能：
再投一次：
同样是要计算总的分数。
再投一次：

就是这样按照套路计算一遍再比较结果即可。

3.只是知道观测序列，想知道骰子是怎么样的？

这个算法在后面会细讲，Baum-welch算法。

④HMM的公式角度

下面正式开始讲解，上面只是过一个印象。

1.马尔科夫模型

要看隐马尔科夫自然先动马尔科夫是什么。已知现在有N个有序的随机变量，根据贝叶斯他们的联合分布可以写成条件连乘：
p(x_1,x_2,x_3,...,x_N) = p(x1|x_2,x_3,...,x_N)p(x_2,x_3,...,x_N)
p(x_2,x_3,...x_N) = p(x_2|x_3,x_4,...x_N)p(x_3,x_4,...x_N)
再继续拆：p(x_1,x_2,x_3...x_N) = \prod_{i=1}^{n}p(x_i|x_{i-1},...,x_1)
马尔科夫链就是指，序列中的任何一个随机变量在给定它的前一个变量的分布于更早的变量无关。也就是说当前的变量只与前一个变量有关，与更早的变量是没有关系的。用公式表示：
p(x_n|x_{n-1},x_{n-2},...,x_1) = p(x_n|x_{n-1})在这个假设的前提下：p(x_1,x_2,x_3,...x_N) = p(x_1)\prod_{i=1}^{n}p(x_i|x_{i-1})这个式子表示的就是一阶马尔科夫模型：

上图表示的是一阶的马尔科夫模型，一阶的意思是只和前一个variable相关，二阶那自然就是和前两个variable相关了，所以M阶的:p(x_n|x_{n-1},...,x_1) = p(x_n|x_{n-1},..x_{n-M})但是这样带来的问题就是指数爆炸的问题，虽然达到了关联更早变量的能力，但是计算能力很大。每一个变量等于是指数级的计算量非常大。那么有没有一种方法可以使得当前变量和更早的变量关联起来呢？又不需要这么多参数？

2.隐马尔科夫模型

于是我们引入了隐变量，做到了一个变量可以和更早的变量关联起来。使用隐变量构成一阶马尔科夫模型，可观测变量和隐变量关联，这样就可以得到一类模型，也就是状态空间模型。

Z_n就是隐变量，x_n为可观测变量。那么x_n,z_n的联合概率：
p(x_1,...,x_n,z_1,...,z_n) = p(x_1|x_2,x_3,...x_n,z_1,z_2,...z_n)p(x_2,x_3,...x_n,z_1,...z_n)
=p(x_1|z_1)...p(x_n,|z_n)p(z_1,...z_n)
=p(z_1)[\prod_{i=1}^{n}p(x_i|z_i)][\prod_{i=2}^{n}p(z_i|z_{i-1})]
于是就被分解成了三个部分：p(z_1),p(z_n|z_{n-1}),p(x_n|z_n)，这三个东西分别是，初始概率，转移概率和发射概率。这个时候，每一个变量之间已经不存在有马尔科夫性了，因为每一个变量x都会通过对应的隐变量和之前的变量相关，所以无法再从p(x_n|x_{n-1},...X_1)中拿到任何一个变量了。其实隐变量其实就像是一个媒介，关联起所有的向量。当z_n是离散的时候，这个模型就是隐马尔科夫模型了。

3.隐马尔科夫模型的定义

根据上面的描述，已经知道了隐马尔科夫模型是由三个部分构成，初始概率，转移概率，发射概率。那自然就对应着三个参数了：\lambda = \{A,B,\pi \}。π是初始矩阵，也就是一个Nx1的，N就是隐变量的数量；A是状态转移矩阵，A = [α_{ij}]_{NxN}，因为这是代表转移到每一个隐变量的概率，包括自己的。B就是发射概率了。
使用I = [i_1,i_2,i_3,...i_T]来代表隐含状态序列，使用O = [o_1,o_2,o_3,...o_T]
隐状态的种类：Q = (q_1,q_2,...q_N)，有N种，每一种分别是q_i。
所有可能观测到的集合：V = (v_1,v_2,...v_M)
最后就是隐马尔科夫模型的两个条件了，这两个条件前面可是提过的，这两个条件在后面大有作为：
①齐次假设：p(i_{t}|i_{t-1},o_{t-1},i_{t-2},o_{t-2}...i_1,o_1) = p(i_t|i_{t-1})
②观测独立性假设：p(o_t|i_t,o_t,i_{t-1},o_{t-1}...i_1,o_1) = p(o_t|i_t)

4.三个问题的公式求解

问题1：给定模型的参数\lambda = ({A, B, \pi})和观测序列O = \{o_1,o_2,o_3,...o_T\}，计算在λ参数下的O观测序列出现的概率是多少？

这个问题，暴力求解，前向后向算法。首先要介绍的就是暴力求解。

暴力求解

这个算法没有什么卵用，只是用于理解这个模型而已。
按照概率公式，列举所有可能的长度为T的状态序列I = \{i_1,i_2,...i_T\}，求各个状态序列I和观测序列O = \{o_1,o_2...,o_T\}的联合概率p(O,I|\lambda)。
还是常规操作，分解了再说：
P(O|\lambda) = \sum_IP(O,I|\lambda) = \sum_IP(O|I,\lambda)P(I|\lambda)
那么现在就是要求解P(I|\lambda),P(O|I,\lambda)。
P(I|\lambda) = P(I_1,I_2,...I_T|\lambda)
=P(I_1|I_2,I_3...I_T,\lambda)P(I_2,I_3...I_T|\lambda)
=P(I_1|\lambda)P(I_2|I_1,\lambda)P(I_3|I_2,\lambda)...P(I_T|I_{T-1},\lambda)
=P(I_1|\lambda)\prod_{i=2}^{T}P(I_i|I_{i-1},\lambda)
=\pi_{i_1}α_{i_1i_2}α_{i_2i_3}...α_{i_{T-1}i_{T}}
这上面用到了齐次假设。
P(O|I,\lambda) = P(O_1,O_2.O_3,...O_T|I_1,I_2,...I_T,\lambda)
=P(O_1|O_2,O_3,...O_T,I_1,I_2,...I_T,\lambda)P(O_2,O_3,...O_T|I_1,I_2,...I_T,\lambda)
=P(O_1|I_1,\lambda)P(O_2,O_3,...O_T|I_1,I_2,...I_T,\lambda)
=\prod_{i=1}^{T}P(O_i|I_i,\lambda)
=b_{i_1o_1}b_{i_2o_2}...b_{i_To_T}
这里就用到了观测独立假设。
所以，P(O|\lambda) = \sum_{i_1,i_2,...i_T}\pi_{i_1}b_{i_1o_1}a_{i_1i_2}b_{i_2o_2}...a_{i_{T-1}i_T}b_{i_To_T}
所以这个就是暴力求解的过程了，没有什么太抽象的东西，都是直接公式的推导代入即可。如果马尔科夫链很长，一般都是做不完的，所以这个算法也就是助于理解而已。

前向算法

在给定时刻，求出现了观测序列O = (O_1,O_2,...O_T)。假设α_{t}(i) = P(O_1,O_2,...O_t,i_t=q_i|\lambda)表示的就是在观测到O=(O_1,O_2,...O_t)同时最后一个状态是q_i的时候的概率。
求初值：α_1(i) = \pi_ib_{io_1}
对于t = 1,2,3...T：α_{t+1}(i)=[\sum_{j=1}^{N}α_t(j)a_{ji}]b_{io_{t+1}}。这个其实蛮好理解的，等于就是：(上一个节点转移的i个状态的和)*i状态到下一个状态的发送概率
最后：P(O|\lambda) = \sum_{i=1}^{N}α_T(i)
前向算法很直观，比较好理解。

后向算法

后向算法，顾名思义就是往前推的了。所以使用β来表示。
β_t(i) = P(O_{t+1},O_{t+2},O_{t+3},...O_T|i_t = q_i,\lambda)
初值：\beta_T(i) = 1
对于t=T-1,T-2,...1：\beta_t(i) = \sum_{j=1}^{N}(a_{ij}b_{jo_{t+1}}\beta_{t+1}(j))
最终：P(O|\lambda) = \sum_{i=1}^N\pi_ib_{io_{1}}\beta_1(i)
直观理解：这个东西有点骚，你得倒着看。首先既然是得到了未来的一个概率，那么要求之前的，肯定要回退到隐状态啊。怎么回退呢？怎么来的就怎么退，隐状态到观测序列用的是B矩阵发射概率，那自然是用B矩阵里面的来回退了，于是乘上b_{jo_{t+1}}回退到状态转移矩阵，然后再用A矩阵回到的上一个节点。但是这样只是做一条路，总共可不止这么多条，于是再加和即可。

这样就应该很直观了。
公式理解：
\beta_t(i) = P(O_{t+1},O_{t+2},O_{t+3},...O_{T}|i_t=q_i,\lambda)
=\sum_jP(O_{t+1},O_{t+2},O_{t+3},..O_T,i_{t+1} = q_j|i_t=q_i,\lambda)
=\sum_jP(O_{t+1},O_{t+2},...O_T|i_{t+1}=q_j,i_t=q_i,\lambda)P(i_{t+1}=q_j|i_t=q_i,\lambda)
=\sum_jP(O_{t+1},O_{t+2},...O_T|i_{t+1}=q_j,\lambda)a_{ij}
=\sum_jP(O_{t+1}|O_{t+2},O_{t+3},...O_T,i_{t+1}=q_j,\lambda)P(O_{t+2},O_{t+3},...O_T|i_{t+1}=q_j,\lambda)a_{ij}
=\sum_jP(O_{t+1}|i_{t+1}=q_j,\lambda)P(O_{t+2},O_{t+3},...O_T|i_{t+1}=q_j,\lambda)a_{ij}
=\sum_jb_{jo_{t+1}}\beta_{t+1}(j)a_{ij}
=\sum_{j=1}^{N}a_{ij}b_{jo_{t+1}}\beta_{t+1}(j)
这个更直观，公式一推就都出来了。这个推导忘了是再哪看到的了。当时我和我的小伙伴们都惊呆了，所以记得比较清楚。

额外补充一点，对于前向概率和后向概率之间的关系：

前向概率和后向概率的关系

我们定义一个概率表达P(O,i_t=q_i|\lambda)，在参数给定了的情况下，求解出现了观测序列O且当前状态是q_i的概率是多少。分解并推导一下这个公式：
P(O,i_t=q_i|\lambda) = P(O|i_t=q_i,\lambda)P(i_t=q_i|\lambda)
=P(O_1,O_2...O_T|i_t=q_i,\lambda)P(i_t=q_i|\lambda)
=P(O_1,...O_t,O_{t+1},...O_T|i_t=q_i,\lambda)P(i_t=q_i|\lambda)
=P(O_1,...O_t|i_t=q_i,\lambda)P(O_{t+1},...O_T|O_1,...O_t,i_t=q_i.\lambda)P(i_t=q_i|\lambda)
=P(O_1,...O_t,i_t=q_i|\lambda)P(O_{t+1},...O_T|i_t=q_i,\lambda)
=α_t(i)\beta_t(i)
记：\gamma_t(i) = P(i_t=q_i|O,\lambda)
\gamma_t(i)=P(i_t=q_i|O,\lambda)
=\frac{P(i_t=q_i,O|\lambda)}{P(O|\lambda)}
=\frac{\alpha_t(i)\beta_t(i)}{\sum_{i=1}^{N}\alpha_t(i)\beta_t(i)}
\gamma_t(i)表示的就是在给定了观测序列和参数的情况下，当前状态是q_i的概率是多少。
还有另外一个比较重要的：给定了观测值和参数，求在t时刻状态是q_i，t+1时刻状态是q_j的概率，记为：ξ_t(i,j) = P(i_t=q_i,i_{t+1}=q_j|O,\lambda)
推导一下，这个在后面baum-welch算法大有作为。
ξ_t(i,j) = P(i_t=q_i,i_{t+1}=q_j|O,\lambda)
=P(i_t=q_i,i_{t+1}=q_j,O|\lambda)/P(O|\lambda)
=P(i_t=q_i,i_{t+1}=q_j,O|\lambda)/\sum_{i=1}^{N}\sum_{j=1}^{N}P(i_t=q_i,i_{t+1}=q_j,O|\lambda)
主要就是P(i_t=q_i,i_{t+1}=q_j,O|\lambda)怎么求。这个和前向算法有点像，直接画图理解：

对比图片其实很明显了：P(i_t=q_i,i_{t+1}=q_j,O|\lambda) = \alpha_t(i)a_{ij}b_{jO_{t+1}}\beta_{t+1}(j)
所以ξ_t(i,j) = \frac{\alpha_t(i)a_{ij}b_{jO_{t+1}}\beta_{t+1}(j)}{\sum_{i=1}^{N}\sum_{j=1}^{N}\alpha_t(i)a_{ij}b_{jO_{t+1}}\beta_{t+1}(j)}
期望：
在观测状态O已知，参数λ已知的情况下，状态i出现的期望：\sum_{t=1}^{T}\gamma_t(i)
在观测状态O已知，参数\lambda已知的情况下，状态i转移到状态j的期望：\sum_{t=1}^{T-1}ξ_t(i,j)

这样的话那么第一个问题就解决了。

问题2：已知观测序列O=(O_1,O_2,...O_T)，估计模型的参数\lambda=(A,B,\pi)的参数，使得在该模型下的P(O|\lambda)最大。

baum-welch算法求解-EM思想

求极大，自然就是EM了。这里直接就使用Q函数：
Q(θ,θ^{(i)}) = \sum_{I}P(I|O,θ^{(i)})logP(O|I,θ)P(I|θ)
还是再提一下，这里直接使用Q函数，这是李航老师里面的正式理解，和通俗理解就差log里面除一个隐变量的分布了。虽然式子不一样，但是求导之后下面这个常数Q分布是可以去掉的，没有任何影响，因为参数使用的是前一次迭代求出来的，所以两者没有区别。
稍微化简一下：
Q(\lambda,\hat\lambda) = \sum_I(lnP(O,I|\lambda))P(I|O,\hat\lambda)
=\sum_I(InP(O,I|\lambda))\frac{P(O,I|\hat\lambda)}{P(O|\hat\lambda)}
由于O = (O_1,O_2,...O_T)是已知的，所以P(O|\hat\lambda)是已知的，那么\sum_I(InP(O,I|\lambda))\frac{P(O,I|\hat\lambda)}{P(O|\hat\lambda)}\propto\sum_I(InP(O,I|\lambda))P(O,I|\hat\lambda)
虽然P(O,I|\hat\lambda)是可以求出来的，但是它是在加和里面，不可以提出来，所以是不能去掉的，而P(O|\hat\lambda)是和加和没有关系的，所以可以正比去掉，反正对整个函数的趋势没有影响的。
前面推导过：P(O,I|\lambda) = \pi_{i_1}b_{i_1o_1}a_{i_1i_2}b_{i_2o_2}...a_{i_{T-1}i_T}b_{i_To_T}
代进去：Q(\lambda,\hat\lambda) = \sum_Iln\pi_{i_1}P(O,I|\hat\lambda) + \sum_I(\sum_{t=1}^{T-1}lna_{i_ti_{t+1}})P(O,I|\hat\lambda)+\sum_I(\sum_{t=1}^{T}lnb_{i_to_t})P(O,I|\hat\lambda)
为什么要分成三个部分呢？因为求导一下就没一边了，开心。
求\pi：
\pi有一个条件：\sum_{i=1}^{N}\pi_i = 1。既然是有条件了，拉格朗日乘子法就用上了。化简一下上式：\sum_Iln\pi_{i_1}P(O,I|\hat\lambda) = \sum_Iln\pi_{i_1}P(O,i_1=i|\hat\lambda)拉格朗日乘子法：
[\sum_{i=1}^{N}ln\pi_{i_1}P(O,i_1=i|\hat\lambda)]+\beta(\sum_{i=1}^{N}\pi_i-1)
求导：
P(O,i_1=1|\hat\lambda) + \beta\pi_1 = 0
P(O,i_1=2|\hat\lambda) + \beta\pi_2 = 0
P(O,i_1=3|\hat\lambda) + \beta\pi_3 = 0
P(O,i_1=4|\hat\lambda) + \beta\pi_4 = 0...
上面那几条式子全部加起来：P(O|\hat\lambda) = -\beta
所以\pi_i = \frac{P(O,i_1=i|\hat\lambda)}{\sum_{i=1}^{N}P(O,i_1=i|\hat\lambda)}=\gamma_1(i)
求A：
这里的条件就是\sum_{j=1}^{N}a_{ij} = 1，还是延用上面的方法，拉格朗日求导，其实就是改一下上面的公式就好了。所以有：
a_{ij} = \frac{\sum_{t=1}^{T-1}P(O,i_t=i,i_{t+1}=j|\hat\lambda)}{\sum_{i=1}^{T-1}P(O,i_t=i|\hat\lambda)}=\frac{\sum_{t=1}^{T-1}ξ_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)}
求B：
这里的条件就是\sum_{j=1}^{M}b_{io_{j}}=1，于是：
b_{ij} = \frac{\sum_{t=1,o_t=v_j}^{T}\gamma_t(i)}{\sum_{t=1}^{T}\gamma_t(i)}
这样就是整个更新算法了。

有监督学习

如果已经有了一堆标注好的数据，那就没有必要再去玩baum-welch算法了。这个算法就很简单了，多的不说，直接上图：

问题3：如果有了参数\lambda和观测序列，求P(I|O,\lambda)的概率最大的隐含序列I是什么。

上面骰子的例子已经举过了，再讲一个下雨天的例子：

day	rain	sun
rain	0.7	0.3
sun	0.4	0.6

转移概率如上

day	walk	shop	clean
rain	0.1	0.4	0.5
sun	0.6	0.3	0.1

发射概率入上
初始概率\pi = (0.6,0.4)
已知模型参数\lambda和观测三天行为(walk,shop,clean)，求天气。
①首先是初始化：
ξ_1 = \pi_{rain}*P(rain|walk) = 0.06
ξ_2 = \pi_{sun}*P(sun|walk) = 0.24
②第二天：
ξ_1 = max[0.06*P(rain|rain),0.024*P(sun|rain)]*P(rain|shop) = 0.0384
φ_1 = 2
ξ_2 = max[0.06*P(rain|sun),0.024*P(sun|sun)]*P(sun|shop) = 0.0432
φ_2 = 2
③最后一天：
ξ_1 = max[0.0384*P(rain|rain),0.0432*P(sun|rain)]*P(rain|clean) = 0.01344
φ_1 = 1
ξ_2 = max[0.0384*P(rain|sun),0.0432*P(sun|sun)]*P(sun|clean) = 0.002592
φ_2 = 2
接着就是回溯了，查看最后一天哪个最大，倒着找，最后就是2,1,1了，也就是(sun,rain,rain)。

⑤Hidden Markov Model代码实现

1.工具类的实现

class Tool(object):
    infinite = float(-2**31)
    @staticmethod
    def log_normalize(a):
        s = 0
        for x in a:
            s += x
        if s == 0:
            print('Normalize error,value equal zero')
            return
        s = np.log(s)
        for i in range(len(a)):
            if a[i] == 0:
                a[i] = Tool.infinite
            else:
                a[i] = np.log(a[i]) - s
        return a

    @staticmethod
    def log_sum(a):
        if not a:
            return Tool.infinite
        m = max(a)
        s = 0
        for t in a:
            s += np.exp(t - m)
        return m + np.log(s)


    @staticmethod
    def saveParameter(pi,A, B, catalog):
        np.savetxt(catalog + '/' + 'pi.txt', pi)
        np.savetxt(catalog + '/' + 'A.txt', A)
        np.savetxt(catalog + '/' + 'B.txt', B)

准备了几个静态方法，一个就是log标准化，也就是把所有的数组都归一化操作，变成一个概率，log算加和，保存矩阵。所有的操作都是使用log表示，如果单单是只用数值表示的话，因为涉及到很多的概率相乘，很多时候就变很小，根本检测不出。而用log之后下限会大很多。

2.有监督算法的实现

其实就是根据上面的公式统计即可。

def supervised(filename):
    '''
    The number of types of lastest variable is four,0B(begin)|1M(meddle)|2E(end)|3S(sigle)
    :param filename:learning fron this file
    :return: pi A B matrix
    '''
    pi = [0]*4
    a = [[0] * 4 for x in range(4)]
    b = [[0] * 65535 for x in range(4)]
    f = open(filename,encoding='UTF-8')
    data = f.read()[3:]
    f.close()
    tokens = data.split('  ')

    #start training
    last_q = 2
    old_process = 0
    allToken = len(tokens)
    print('schedule : ')
    for k, token in enumerate(tokens):
        process = float(k) /float(allToken)
        if process > old_process + 0.1:
            print('%.3f%%' % (process * 100))
            old_process = process
        token = token.strip()
        n = len(token)
        #empty we will choose another
        if n <= 0:
            continue
        #if just only one
        if n == 1:
            pi[3] += 1
            a[last_q][3] += 1
            b[3][ord(token[0])] += 1
            last_q = 3
            continue
        #if not
        pi[0] += 1
        pi[2] += 1
        pi[1] += (n-2)
        #transfer matrix
        a[last_q][0] += 1
        last_q = 2
        if n == 2:
            a[0][2] += 1
        else:
            a[0][1] += 1
            a[1][1] += (n-3)
            a[1][2] += 1
        #launch matrix
        b[0][ord(token[0])] += 1
        b[2][ord(token[n-1])] += 1
        for i in range(1, n-1):
            b[1][ord(token[i])] += 1
    pi = Tool.log_normalize(pi)
    for i in range(4):
        a[i] = Tool.log_normalize(a[i])
        b[i] = Tool.log_normalize(b[i])
    return pi, a, b

按照公式来，一个一个单词判断，如果是单个的那自然就是single，所以当n == 1的时候直接就是\pi[3] += 1。初始状态last_q = 2是因为一开始肯定是结束之后才接着开始的，所以自然就是2，end。之后都是比较常规了。

3.viterbi算法的实现

def viterbi(pi, A, B, o):
    '''
    viterbi algorithm
    :param pi:initial matrix
    :param A:transfer matrox
    :param B:launch matrix
    :param o:observation sequence
    :return:I
    '''
    T = len(o)
    delta = [[0 for i in range(4)] for t in range(T)]
    pre = [[0 for i in range(4)] for t in range(T)]
    for i in range(4):
        #first iteration
        delta[0][i] = pi[i] + B[i][ord(o[0])]
    for t in range(1, T):
        for i in range(4):
            delta[t][i] = delta[t-1][0] + A[0][i]
            for j in range(1, 4):
                vj = delta[t-1][j] + A[0][j]
                if delta[t][i] < vj:
                    delta[t][i] = vj
                    pre[t][i] = j
            delta[t][i] += B[i][ord(o[t])]
    decode = [-1 for t in range(T)]
    q = 0
    for i in range(1, 4):
        if delta[T-1][i] > delta[T-1][q]:
            q = i
    decode[T-1] = q
    for t in range(T-2, -1, -1):
        q = pre[t+1][q]
        decode[t] = q
    return decode

根据上面的例子来就好，先找到转移概率最大的，迭代到最后找到概率最大的序列之后倒着回来找即可。最后就得到一串编码，然后使用这段编码来进行划分。

def segment(sentence, decode):
    N = len(sentence)
    i = 0
    while i < N:
        if decode[i] == 0 or decode[i] == 1:
            j = i+1
            while j < N:
                if decode[j] == 2:
                    break
                j += 1
            print(sentence[i:j+1],"|",end=' ')
            i = j+1
        elif decode[i] == 3 or decode[i] == 2:  # single
            print (sentence[i:i + 1], "|", end=' ')
            i += 1
        else:
            print ('Error:', i, decode[i] , end=' ')
            i += 1

这个就是根据编码划分句子的函数了。首先是通过有监督的学习得到参数\lambda = (A, B, \pi)，然后用viterbi算法得到编码序列，再用segment函数做划分即可。

if __name__ == '__main__':
    pi = np.loadtxt('supervisedParam/pi.txt')
    A = np.loadtxt('supervisedParam/A.txt')
    B = np.loadtxt('supervisedParam/B.txt')
    f = open("../Data/novel.txt" , encoding='UTF-8')
    data = f.read()[3:]
    f.close()
    decode = viterbi(pi, A, B, data)
    segment(data, decode)

执行代码。

效果当然可以了，毕竟是有监督。无监督的就没有这么秀了。

4.baum-welch算法的实现

参考上面三个公式，除了B有点难更新之外其他的都很简单。

def cal_alpha(pi, A, B, o, alpha):
    print('start calculating alpha......')
    for i in range(4):
        alpha[0][i] = pi[i] + B[i][ord(o[0])]
    T = len(o)
    temp = [0 for i in range(4)]
    del i
    for t in range(1, T):
        for i in range(4):
            for j in range(4):
                temp[j] = (alpha[t-1][j] + A[j][i])
            alpha[t][i] = Tool.log_sum(temp)
            alpha[t][i] += B[i][ord(o[t])]
    print('The calculation of alpha have been finished......')

def cal_beta(pi, A, B, o, beta):
    print('start calculating beta......')
    T = len(o)
    for i in range(4):
        beta[T-1][i] = 1
    temp = [0 for i in range(4)]
    del i
    for t in range(T-2, -1, -1):
        for i in range(4):
            beta[t][i] = 0
            for j in range(4):
                temp[j] = A[i][j] + B[j][ord(o[t + 1])] + beta[t + 1][j]
            beta[t][i] += Tool.log_sum(temp)
    print('The calculation of beta have been finished......')

def cal_gamma(alpha, beta, gamma):
    print('start calculating gamma......')
    for t in range(len(alpha)):
        for i in range(4):
            gamma[t][i] = alpha[t][i] + beta[t][i]
        s = Tool.log_sum(gamma[t])
        for i in range(4):
            gamma[t][i] -= s
    print('The calculation of gamma have been finished......')

def cal_kesi(alpha, beta, A, B, o, ksi):
    print('start calculating ksi......')
    T = len(o)
    temp = [0 for i in range(16)]
    for t in range(T - 1):
        k = 0
        for i in range(4):
            for j in range(4):
                ksi[t][i][j] = alpha[t][i] + A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j]
                temp[k] = ksi[t][i][j]
                k += 1
        s = Tool.log_sum(temp)
        for i in range(4):
            for j in range(4):
                ksi[t][i][j] -= s
    print('The calculation of kesi have been finished......')

α,β,ξ,γ的更新按照公式即可。

def update(pi, A, B, alpha, beta, gamma, ksi, o):
    print('start updating......')
    T = len(o)
    for i in range(4):
        pi[i] = gamma[0][i]
    s1 = [0 for x in range(T-1)]
    s2 = [0 for x in range(T-1)]
    for i in range(4):
        for j in range(4):
            for t in range(T-1):
                s1[t] = ksi[t][i][j]
                s2[t] = gamma[t][i]
            A[i][j] = Tool.log_sum(s1) - Tool.log_sum(s2)
    s1 = [0 for x in range(T)]
    s2 = [0 for x in range(T)]
    for i in range(4):
        for k in range(65536):
            if k%5000 == 0:
                print(i, k)
            valid = 0
            for t in range(T):
                if ord(o[t]) == k:
                    s1[valid] = gamma[t][i]
                    valid += 1
                s2[t] = gamma[t][i]
            if valid == 0:
                B[i][k] = -Tool.log_sum(s2)
            else:
                B[i][k] = Tool.log_sum(s1[:valid]) - Tool.log_sum(s2)
    print('baum-welch algorithm have been finished......')

最后再封装一把：

def baum_welch(pi, A, B, filename):
    f = open(filename , encoding='UTF-8')
    sentences = f.read()[3:]
    f.close()
    T = len(sentences)   # 观测序列
    alpha = [[0 for i in range(4)] for t in range(T)]
    beta = [[0 for i in range(4)] for t in range(T)]
    gamma = [[0 for i in range(4)] for t in range(T)]
    ksi = [[[0 for j in range(4)] for i in range(4)] for t in range(T-1)]
    for time in range(1000):
        print('Time : ', time)
        sentence = sentences
        cal_alpha(pi, A, B, sentence, alpha)
        cal_beta(pi, A, B, sentence, beta)
        cal_gamma(alpha, beta, gamma)
        cal_kesi(alpha, beta, A, B, sentence, ksi)
        update(pi, A, B, alpha, beta, gamma, ksi, sentence)
        Tool.saveParameter(pi, A, B, 'unsupervisedParam')
        print('Save matrix successfully!')

初始化矩阵：

def inite():
    pi = [random.random() for x in range(4)]
    Tool.log_normalize(pi)
    A = [[random.random() for y in range(4)] for x in range(4)]
    A[0][0] = A[0][3] = A[1][0] = A[1][3] = A[2][1] = A[2][2] = A[3][1] = A[3][2] = 0
    B = [[random.random() for y in range(65536)] for x in range(4)]
    for i in range(4):
        A[i] = Tool.log_normalize(A[i])
        B[i] = Tool.log_normalize(B[i])
    return pi , A , B

最后跑一遍就OK了。
无监督算法使用的是EM框架，肯定会存在局部最大的问题和初始值敏感，跑了56次，用的谷歌GPU跑，代码没有经过优化，跑的贼慢。

最后的效果也是惨不忍睹。

总结

概率图模型大多都是围绕着三个问题展开的，求观测序列的概率最大，求隐含序列的概率最大，求参数。MEMM,RCF大多都会围绕这几个问题展开。求观测序列的概率，暴力求解是为了理解模型，前向后向算法才是真正有用的；概率最大的隐含序列viterbi算法，动态规划的思想最后回溯回来查找；求参数就是套用EM框架求解。
HMM是属于生成模型，首先求出P(O|\lambda)，转移概率和表现概率直接建模，也就是对样本密度建模，然后才进行推理预测。事实上有时候很多NLP问题是和前后相关，而不是只是和前面的相关，HMM这里明显是只和前面的隐变量有关，所以还是存在局限性的。对于优化和模型优缺点等写了MEMM和RCF再一起总结。

最后附上GitHub代码：

https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/HMM