机器学习（十）——概率图模型之隐马尔可夫模型

作者: 夏普123 | 来源:发表于2021-12-12 15:38 被阅读0次

嗯，终于迎来了终极大Boss。

1. 概论

想要理解清楚隐马尔可夫模型确实比之前要难一些，但是我尽量讲清楚。先来看一个问题：
假设天气的状况分为：晴天、多云、雨天。
我想预测明天的天气状况。
这里我们假设：每天的天气状况只跟前几天的天气状况有关系，而不去考虑其他影响因素，例如风力、气压等等。—— 这就是著名的马尔可夫假设。即：假设模型的当前状态仅仅依赖于前面的几个状态。——当然，这种假设非常粗糙，并且因此可能将一些非常重要的信息丢失，但是它缺极大的简化了问题。
一个马尔科夫过程是状态间的转移仅依赖于前n个状态的过程。这个过程被称之为n阶马尔科夫模型，其中n是影响下一个状态选择的（前）n个状态。最简单的马尔科夫过程是一阶模型，它的状态选择仅与前一个状态有关。为了把问题更加简化，我们后面都将基于一阶马尔可夫过程来讨论。

根据一阶马尔可夫过程，明天的天气状况，只依赖于今天的天气状况。由于从晴天到多云，从晴天到雨天，从晴天到晴天的变化不是确定的，也就是说，晴天的后面不一定还是晴天，也不一定是雨天，也不一定是多云，任何一个状态都有可能是所有状态的下一个转移状态，我们假设每一个状态转移都有一个概率值，称为状态转移概率——这是从一个状态转移到另一个状态的概率，那么我们有多少种状态转移的可能性呢？按上面的例子，我们总共有如下几种：
晴天->晴天；
晴天->雨天；
晴天->多云；
雨天->晴天；
雨天->雨天；
雨天->多云；
多云->晴天；
多云->雨天；
多云->多云；
所以，总共是 $M*M=M^2$ 种状态转移过程。我们把上面这几种过程的概率都列出来，就形成了一个状态转移概率矩阵，如下所示：

p	晴天	多云	雨天
晴天	0.5	0.375	0.125
多云	0.25	0.125	0.625
雨天	0.25	0.375	0.375

注意，每一行的总和为1。第一行表示，今天是晴天，那么明天也是晴天的概率为0.5，今天是晴天，明天是多云的概率是0.375，今天是晴天，明天是雨天的概率是0.125。
上面的概率转移矩阵，我们用数学公式来表示为：
$A = \begin{bmatrix} 0.5 & 0.25 & 0.25 \\ 0.375 & 0.125 & 0.375 \\ 0.125 & 0.625 & 0.375 \end{bmatrix}$

注意，此时是按列来看的。第一列代表今天是晴天，明天是晴天、多云、雨天的概率，依次类推。

假设今天是晴天，我们初始化一个向量，来表示今天是晴天：
$pi=\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$ ，

那么
$A*pi = \begin{bmatrix} 0.5 & 0.25 & 0.25 \\ 0.375 & 0.125 & 0.375 \\ 0.125 & 0.625 & 0.375 \end{bmatrix} \times \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} (0.5*1+0.25*0+0.25*0) \\(0.375*1+0.125*0+0.375*0)\\(0.125*1+0.625*0+0.375*0)\end{bmatrix}=\begin{bmatrix} 0.5 \\ 0.375 \\ 0.125 \end{bmatrix}$
就代表明天是晴天、多云、雨天的概率。

有了明天的天气概率，我们要预测后天的概率要怎么算呢？显然，后天的概率就是用明天的天气概率去乘以概率转移矩阵，也即：
$后天的概率=明天的概率*A=A*pi*A=A^2*pi$
同理，大后天的天气概率为：
$大后天的天气概率=后天的天气概率*A=A^3*pi$

可以看到，我们如果要预测未来n天的天气，我们是通过今天的天气一步一步向前预测推导的。后面每一天的天气状态都只跟前一天的天气状态相关。如此形成了一个链条，这就是马尔可夫链。

好了，现在我们再来看什么是隐马尔可夫链。

假设有一个被关在地牢里的家伙，他无法看到外面的天气状况。但是他可以看到地牢里的干湿状态：干燥、稍干、潮湿、湿润。地牢的干湿状态与天气的状况有一定的概率关系。那么这个被关在地牢里的家伙他能否预测未来的天气情况呢？
我们把地牢的干湿程度记为观测状态（用O表示），对于地牢的人来说，天气状态就是隐藏状态（用S表示），用图表示为：

这就是隐马尔可夫模型。

从上面的介绍中，我们知道，一个隐马尔可夫模型由如下五个元素：

隐含状态 S。这些状态之间满足马尔可夫性质，是马尔可夫模型中实际所隐含的状态。这些状态通常无法通过直接观测而得到。
（例如上面的天气状态）
可观测状态 O。在模型中与隐含状态相关联，可通过直接观测而得到。（例如上面的地牢状态）
初始状态概率矩阵 $\pi$ 。表示隐含状态在初始时刻t=1的概率矩阵
隐含状态转移概率矩阵 A。
观测状态转移概率矩阵 B。有的地方也叫混淆矩阵，表示给定一个隐藏状态后得到的观察状态的概率。

一般的，可以用 $\lambda=(A,B,\pi)$ 三元组来简洁的表示一个隐马尔可夫模型。隐马尔可夫模型实际上是标准马尔可夫模型的扩展，添加了可观测状态集合和这些状态与隐含状态之间的概率关系。

对于隐马尔可夫模型，我们根据实际情况的不同，一般会遇到3类问题：
1) 评估模型
给定观测序列 $O=O_1,O_2,O_3,…,O_t$ 和模型参数 $\lambda=(A,B,\pi)$ ，怎样有效计算某一观测序列的概率，进而可对该HMM做出相关评估。例如，已有一些模型参数各异的HMM，给定观测序列 $O=O_1,O_2,O_3,…,O_t$ ，我们想知道哪个HMM模型最可能生成该观测序列。通常我们利用forward算法分别计算每个HMM产生给定观测序列O的概率，然后从中选出最优的HMM模型。
这类评估的问题的一个经典例子是语音识别。在描述语言识别的隐马尔科夫模型中，每个单词生成一个对应的HMM，每个观测序列由一个单词的语音构成，单词的识别是通过评估进而选出最有可能产生观测序列所代表的读音的HMM而实现的。
2) 解码问题
给定观测序列 $O=O_1,O_2,O_3,…,O_t$ 和模型参数 $\lambda=(A,B,\pi)$ ，怎样寻找某种意义上最优的隐状态序列。在这类问题中，我们感兴趣的是马尔科夫模型中隐含状态，这些状态不能直接观测但却更具有价值，通常利用Viterbi算法来寻找。
这类问题的一个实际例子是中文分词，即把一个句子如何划分其构成才合适。例如，句子“发展中国家”是划分成“发展-中-国家”，还是“发展-中国-家”。这个问题可以用隐马尔科夫模型来解决。句子的分词方法可以看成是隐含状态，而句子则可以看成是给定的可观测状态，从而通过建HMM来寻找出最可能正确的分词方法。
3) 学习问题
即HMM的模型参数 $\lambda=(A,B,\pi)$ 未知，如何调整这些参数以使观测序列 $O=O_1,O_2,O_3,…,O_t$ 的概率尽可能的大。通常使用Baum-Welch算法以及Reversed Viterbi算法解决。

2. 评估模型

我们先来看第一个问题，即评估模型。评估模型主要解决的问题是：给定模型 $\lambda=(A,B,\pi)$ 和观测序列 $O={O_1,O_2,O_3,…,O_t}$ ，计算模型 $\lambda$ 下观测到序列O出现的概率 $P(O|\lambda)$ 。
对于求解这个问题，我们可以使用穷举法。还是以上面那个例子，天气的状态转移矩阵我们已经知道了，为：
$A = \begin{bmatrix} 0.5 & 0.25 & 0.25 \\ 0.375 & 0.125 & 0.375 \\ 0.125 & 0.625 & 0.375 \end{bmatrix}$
假设，这三种天气导致地牢干燥和潮湿的概率为（为了简化问题，我们这里只假设地牢只有两种状态：干燥和潮湿）：
$B = \begin{bmatrix} 0.8 & 0.6 & 0.1 \\ 0.2 & 0.4 & 0.9 \end{bmatrix}$
如果我们要求地牢的干湿程度为O={干燥、潮湿、潮湿}的概率，我们要这么计算：
如果三天都为晴天，得到上述观测序列的概率为 $P_1=P(干燥|晴天）*P(潮湿|晴天)*P(潮湿|晴天)$
如果三天的天气状况为：晴天，晴天，多云，则得到上面观测序列的概率为： $P_2=P(干燥|晴天）*P(潮湿|晴天)*P(潮湿|多云)$
如果三天的天气状况为：晴天，晴天，雨天，则得到上面观测序列的概率为： $P_3=P(干燥|晴天）*P(潮湿|晴天)*P(潮湿|雨天)$
如果三天的天气状况为：晴天，多云，晴天，则得到上面观测序列的概率为： $P_4=P(干燥|晴天）*P(潮湿|多云)*P(潮湿|晴天)$
……
最后的得到概率应该为上面穷举的所有概率之和。
当然，这里我们还没有计算分别得到上面这所有天气情况下的概率。这样计算的话，将会非常复杂。在观察值较多的情况下，显然是不太现实的。因此我们要用到《数据驱动的决策分析》中学到的动态规划思想。

动态规划
动态规划的思想，其实就是一个递归思想。

如果用穷举法，则从A到E一共有3×3×2=18条不同的路径，逐个计算每条路径的长度，总共需要进行4×18=72次加法计算；对18条路径的长度做两两比较，找出其中最短的一条，总共要进行18-1=17次比较。
如果从A到C的站点有k个，则总共有3k-1×2条路径，用穷举法求最短路径总共要进行(k+1)3k-1×2次加法，3k-1×2-1次比较。当k的值增加时，需要进行的加法和比较的次数将迅速增加。例如当k=10时，加法次数为433,026次，比较39,365次。
显然，当问题的规模很大时，穷举法无法解决这样的问题。
动态规划用一种全新的思路来考虑这样的问题。
首先，把一个复杂的问题归结为若干个与原问题性质完全相同，但规模较小的子问题，每一个子问题又可以递归地归结为性质相同，规模更小的下层子问题。如此不断进行，一直到子问题非常简单，可以解决为止。
然后，从这些最简单的子问题出发，依次计算上一层子问题的最优解，直到求出原问题的最优解。

以最短路径问题为例，来说明动态规划的基本思想。
从A到E的最短路径S(A)，可以转化为三个性质完全相同，但规模较小的子问题，即分别从B1、B2、B3到E的最短路径，记为S(B1), S(B2), S(B3)。

同样，计算S(B1)、 S(B2) 、S(B3)又可以归结为性质完全相同，但规模更小的问题，即分别求C1，C2，C3到E的最短路径问题S(Ci) (i=1, 2, 3)。

计算S(Ci)又可以归结为求从D1和D2到E的最短路径S(D1)和S(D2)这两个子问题。

S(D1)和S(D2)是以知的，它们分别是：S(D1)=5，S(D2)=2。因而，可以从这两个值开始，逆向递归计算S(A)的值。

向前算法
现在回到我们的问题。我们现在要评估观测数据O出现的概率。我们可以这样分解：
先看第一天地牢为干燥的概率。因为我们有三种天气状态，每种天气状态都有可能会导致地牢干燥。因此第一天地牢干燥的概率为：
$P(O_1|\lambda)=\pi_1*b_1(O_1)+\pi_2*b_2(O_1)+\pi_3*b_3(O_1)$
假设： $pi=\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}$ ，也即第一天是晴天，那么 $\pi_1=1,\pi_2=0,\pi_3=0$ 。根据上面的状态转移矩阵，我们知道， $b_1=0.8$ ，代表晴天到干燥的概率：P(干燥|晴天)， $b_2=0.6，b_3=0.1$ 分别代表P(干燥|多云)，P(干燥|雨天)。
我们记： $\alpha_1(1)=\pi_1*b_1(O_1),\alpha_2(1)=\pi_2*b_2(O_1),\alpha_3(1)=\pi_3*b_3(O_1)$
因此， $P(O_1|\lambda)=\alpha_1(1)+\alpha_2(1)+\alpha_3(1)=1*0.8+0*0.6+0*0.6=0.8$

再来看第二天地牢为潮湿的概率。
当 $t=2$ 时，输出为 $O_1O_2$ 。假设第二天是晴天（记为 $s_2=q_1$ ），那么地牢为潮湿的概率为：
$P(O_1O_2,s_2=q_1|\lambda)=\pi_1*b_1(O_1)*a_{11}*b_1(O_2)+\pi_2*b_2(O_1)*a_{21}*b_1(O_2)+\pi_3*b_3(O_1)*a_{31}*b_1(O_2)\\=\alpha_1(1)*a_{11}*b_1(O_2)+\alpha_2(1)*a_{21}*b_1(O_2)+\alpha_3(1)*a_{31}*b_1(O_2)$ ，记为： $\alpha_1(2)$
注意，这里的 $a_{11}、a_{21}、a_{31}$ 是状态转移矩阵A中的元素，代表从晴天、多云、雨天转移到晴天的概率，不是我们上面定义的 $\alpha$ 。
同理，假设第二天是多云，那么地牢为潮湿的概率为：
$P(O_1O_2,s_2=q_2|\lambda)=\alpha_1(1)*a_{12}*b_2(O_2)+\alpha_2(1)*a_{22}*b_2(O_2)+\alpha_3(1)*a_{32}*b_2(O_2)$ ，记为： $\alpha_2(2)$
假设第二天是雨天，那么地牢为潮湿的概率为：
$P(O_1O_2,s_2=q_3|\lambda)=\alpha_1(1)*a_{13}*b_3(O_2)+\alpha_2(1)*a_{23}*b_3(O_2)+\alpha_3(1)*a_{33}*b_3(O_2)$ ，记为： $\alpha_3(2)$
那么，第二天为潮湿的概率就为：
$P(O_1O_2|\lambda)=\alpha_1(2)+\alpha_2(2)+\alpha_3(2)$

同理，我们来看第三天为潮湿的概率
当 $t=3$ 时，输出为 $O_1O_2O_3$ 。
当第三天为晴天时：
$\alpha_1(3)=\alpha_1(2)*a_{11}*b_1(O_3)+\alpha_2(2)*a_{21}*b_1(O_3)+\alpha_3(2)*a_{31}*b_1(O_3)$
当第三天为多云时：
$\alpha_2(3)=\alpha_1(2)*a_{12}*b_2(O_3)+\alpha_2(2)*a_{22}*b_2(O_3)+\alpha_3(2)*a_{32}*b_2(O_3)$
当第三天为多云时：
$\alpha_3(3)=\alpha_1(2)*a_{13}*b_3(O_3)+\alpha_2(2)*a_{23}*b_3(O_3)+\alpha_3(2)*a_{33}*b_3(O_3)$
最后，我们得到最终的概率：
$P(O_1O_2O_3|\lambda)=\alpha_1(3)+\alpha_2(3)+\alpha_3(3)$

由此，我们得到向前算法的一般性步骤（公式）：

step1: 初始化 $\alpha_i(1)=\pi_i*b_i(O_1)$ 。（i表示隐藏状态的计数，例如上面的天气状态：晴天，多云，雨天）

step2: 计算 $\alpha_i(t)=(\sum_{j=1}^N\alpha_j(t-1)*a_{ji})*b_i(O_t)$

step3: $P(O|\lambda)=\sum_{i=1}^N\alpha_i(t)$

3. 解码问题

给定观测序列 $O=O_1,O_2,O_3,…,O_t$ 和模型参数 $\lambda=(A,B,\pi)$ ，找到最可能的状态序列 $I^*=\{i_1^*,i_2^*...i_n^* \}$
同上面的思路，我们可以使用穷举法，列出来所有可能导致地牢干燥、潮湿、潮湿的天气状态组合以及产生地牢干燥、潮湿、潮湿状态的概率。（在上面一节中提到的 $P_1、P_2、P_3、P_4……$ ），然后通过比较这些概率的大小，找到概率最大的那个天气状态组合即可。
显然，当数据量较大时，穷举法在时间上也是不可接受的。同样，我们也考虑使用递归算法来解决上述问题。这就是维比特算法。

维比特算法
维比特算法的思路跟上一节中，我们介绍的动态规划中的思路是一致的，只不过在动态规划中，我们求的是路径最小值，现在我们求的是概率最大值。因此这里就不再赘述了。直接给出步骤和公式：

step1: 初始化。即计算在t=1的时刻，达到某状态的最大概率（这个局部概率我们用 $\delta$ 来表示）。i表示隐藏状态的计数，例如上面的天气状态：晴天，多云，雨天： $\delta_1(i)=\pi_i*b_i(O_1)$ 。

step2: 递推。计算t>1时刻的局部概率。在一阶马尔可夫假设下，状态X在一个状态序列后发生的概率只取决于之前的一个状态。还是以上面的那个例子来说，在t时刻，天气状态为X的概率，只取决于t-1时刻的状态。对于t-1时刻的状态，我们有三种：晴天、多云、雨天。那么，从晴天到达t时刻X的概率为： $P(到达t-1时刻为晴天的最有可能的概率)*P(X|晴天)$
在t时刻X状态下的观察概率为： $P(到达t-1时刻为晴天的最有可能的概率)*P(X|晴天)*P(观察状态|X)$
同样的，从雨天到达t时刻X的观察概率为： $P(到达t-1时刻为雨天的最有可能的概率)*P(X|雨天)*P(观察状态|X)$
所以， $\delta_t(i)=max[\delta_{t-1}(j)a_{ji}]b_i(O_t)$

step3: 计算 $\delta_t(i)$ 的最大值。获得此最大值的观测序列就是我们要找的观测序列。

4. 学习问题

好，终于来到了我们要重点关注，也就是机器学习中经常用到，但是也是最复杂的一个问题，就是已知观测序列 $O=O_1,O_2,O_3,…,O_t$ ，估计模型参数 $\lambda=(A,B,\pi)$ ，使得在该模型参数下出现该观测序列的概率最大。
学习问题要根据观测序列来推导模型参数，这一类的问题对应到概率论中的极大似估计问题。但是这里是有隐含变量的极大似然估计，因此直接无法通过直接求导进行求解，而要通过EM算法来求解这一类问题。
EM算法是一类算法，用于解决有隐含变量的概率模型参数的极大似然估计，具体到隐马尔可夫模型中的具体算法是 Baum-Welch 算法（也叫向前-向后算法）。
直观地理解EM算法，它也可被看作为一个逐次逼近算法：事先并不知道模型的参数，可以随机的选择一套参数或者事先粗略地给定某个初始参数λ0 ，确定出对应于这组参数的最可能的状态，计算每个训练样本的可能结果的概率，在当前的状态下再由样本对参数修正，重新估计参数λ ，并在新的参数下重新确定模型的状态，这样，通过多次的迭代，循环直至某个收敛条件满足为止，就可以使得模型的参数逐渐逼近真实参数。
由于这个算法实在太复杂，这里就不细展开了。感兴趣的同学可以在网上找一些相资料进行学习。

如果是要做一个调包侠，个人觉得只看概论部分就应该可以了。
好了，还是老规矩，最后一部分，让我们用R代码来结束本章内容，也是机器学习系列的最后内容。

5.R代码的实现

例子1: 预测基因序列中的启动因子序列。
这个例子是：首先给定一堆已经区分好了启动子和非启动子的基因序列，我们要找到某种规律或者模型，使得下次我们再提供一种基因序列时，这个模型能够很好的告诉我们，该序列是否属于启动子序列。
这个例子中，不存在隐藏状态。我们为了使用到HMM模型，人为定义了一个隐藏状态。该隐藏状态的序列和观测状态的序列一样，都是样本中的基因序列。且，从隐藏状态到观测状态，相同的碱基之间转换的概率为1，不同的碱基之间转换的概率为0。
另，基因序列是由四个不同的碱基通过不同的排列组合构成。即，一个基因序列有四个状态，分别为：a、c、g、t。为了区分基因序列的起止，我们认为加入了两个状态。S代表基因序列的开始，X代表基因序列的结束。对于初始状态，我们认为只有S。因此初始状态矩阵为c=(1,0,0,0,0,0)。
现在HMM模型中的五个要素，我们已经有了四个。剩下的，我们只需要计算出来隐藏状态转移矩阵A即可。
有了模型之后，我们只需要使用向前算法，即可计算出观测序列出现的概率。理论上，我们如果使用的非启动因子估计出的模型参数，那么在非启动因子的基因序列上，出现该基因序列的概率会比较大，而对于启动因子的基因序列，出现该基因序列的概率会比较小。由此我们可以使用训练数据，采用交叉检验法来衡量模型的好坏。
好了，以上就是下面这个例子的总体思想。现在我们来看代码：

#预测启动子基因序列。
promoters <- read.csv("study/graduate/机器学习/promoters.data", header=F, dec=",", strip.white=TRUE, stringsAsFactors = FALSE)
#V1:+表示启动子，-表示非启动子
#V2:特定序列的识别符
#V3: 核苷酸序列本身（即观测序列，同时也是隐藏序列）
promoters[1,]

#数据拆分成启动子和非启动子。
#subset表示根据条件来筛选子集。这里的3表示选择第3列，丢弃掉1，2列
positive_observations <- subset(promoters, V1=='+', 3)
negative_observations <- subset(promoters, V1=='-', 3)

#在每个序列之前加上S，表示数据的开始，在每个结尾加上X，表示数据的结束
positive_observations<-sapply(positive_observations, function(x) paste("S",x,"X",sep=""))
negative_observations<-sapply(negative_observations, function(x) paste("S",x,"X",sep=""))
positive_observations[1]

#将一串字符串，按每个字母分开,得到一个列表
positive_observations<-strsplit(positive_observations,"")
negative_observations<-strsplit(negative_observations,"")
head(positive_observations[1])

#定义隐藏状态。我们的隐藏状态有以下6个（S,X,a,c,g,t）
states <- c("S", "X", "a", "c", "g", "t")
#定义观测状态。我们的观测状态也有6个
symbols <- c("S", "X", "a", "c", "g", "t")
#定义初始概率，也即pi，表示刚开始的时候，只有S，即从S开始。
startingProbabilities<-c(1,0,0,0,0,0)
#定义发射概率，也即观测状态转移概率矩阵 B。这里是一个对角矩阵，对角线都为1，其余为0
#代表只能从S到S，从a到a，从c到c，从g到g，从t到t
emissionProbabilities<-diag(6)
colnames(emissionProbabilities)<-states #列名为状态值
rownames(emissionProbabilities)<-symbols #行名为观测值
emissionProbabilities

#根据我们的样本，来计算转移概率矩阵。
calculateTransitionProbabilities <- function(data, states) {
  #初始化转移状态矩阵，刚开始都为0
    transitionProbabilities<-matrix(0,length(states),length(states))
    colnames(transitionProbabilities)<-states
    rownames(transitionProbabilities)<-states
  #按顺序循环数据中的每一个字母，统计从状态A状态B的次数，并放到矩阵中。即汇总出每一个状态转移的总数
    for (index in 1:(length(data)-1)) {
        current_state <- data[index]
        next_state <- data[index+1]
        transitionProbabilities[current_state, next_state] <- transitionProbabilities[current_state, next_state] + 1
    }
    #根据次数矩阵，使用sweep函数来计算概率。
    #swepp中的第二个参数1，代表按行计算；2，代表按列计算
    #rowSums函数代表按行进行求和
    #整个函数的意思就是，用每一行的元素除以这一行的数据总和，得到归一化的概率；
    transitionProbabilities<-sweep(transitionProbabilities, 1, rowSums(transitionProbabilities), FUN="/")
    return(transitionProbabilities)
}

#这里相当于把negative_observations从原来的列表转换为了一个向量
#即，把原来每一行数据，都进行了连接，形成了一个长向量，向量长度为原来的行数*每行的字母个数
negative_observation<-Reduce(function(x,y) c(x, y), negative_observations, c())
#使用非启动子基因序列的样本来训练转移概率矩阵。
(transitionProbabilitiesNeg  <- calculateTransitionProbabilities(negative_observation, states))

library("HMM")
#使用HMM的五个参数来初始化隐马尔可夫模型
negative_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesNeg, emissionProbs=emissionProbabilities)

#交叉检验。每次我们从启动子基因序列的样本中拿出一条数据作为检验数据，剩下的作为训练集数据
#先用剩下的数据来构建HMM模型，此时的模型是根据启动子基因序列来构建的
#然后，我们用剩下的这条检验数据，分别带入到之前negative_hmm模型和现在的positive_hmm模型中，采用向前算法来计算出现该基因序列的概率
#理论上，因为这条检验数据是来自于启动子基因序列，因此，通过negative_hmm模型计算出来的概率应该小于使用positive_hmm模型计算出来的概率。
#统计错误率，我们即可知道模型的好坏
incorrect<-0 #定义错误次数统计变量
for (obs in 1:length(positive_observations)) {
  #使用启动子基因序列，构造启动子基因序列模型，方便后面进行对比  
  positive_observation<-Reduce(function(x,y) c(x, y), positive_observations[-obs], c())
  transitionProbabilitiesPos  <- calculateTransitionProbabilities(positive_observation, states)
  positive_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesPos, emissionProbs=emissionProbabilities)

  #获取检验数据
  test_observation<-positive_observations[[obs]]
  final_index<-length(test_observation)
  
  #将检验数据分别带入到两个模型中，使用向前算法，分别计算使用这两个模型得到该序列的概率
  pos_probs<-exp(forward(positive_hmm,test_observation))
  neg_probs<-exp(forward(negative_hmm,test_observation))
  
  #计算出现该序列的概率之和。（上面得到的是序列中每个碱基出现的概率）
  pos_seq_prob<-sum(pos_probs[,final_index])
  neg_seq_prob<-sum(neg_probs[,final_index])
  
  #如果使用positive模型得到的概率小于使用negative模型得到概率，则错误数+1
  if (pos_seq_prob<neg_seq_prob) incorrect<-incorrect+1
}

#反过来，使用positive数据构建模型，使negative数据进行检验。因此以下代码不再做解释
positive_observation <- Reduce(function(x,y) c(x, y), positive_observations, c())
transitionProbabilitiesPos  <- calculateTransitionProbabilities(positive_observation, states)
positive_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesPos, emissionProbs=emissionProbabilities)
  for (obs in 1:length(negative_observations)) {
  negative_observation<-Reduce(function(x,y) c(x, y), negative_observations[-obs], c())
  transitionProbabilitiesNeg <- calculateTransitionProbabilities(negative_observation, states)
  negative_hmm = initHMM(states, symbols, startProbs=startingProbabilities, transProbs=transitionProbabilitiesNeg, emissionProbs=emissionProbabilities)
  
  test_observation<-negative_observations[[obs]]
  final_index<-length(test_observation)
  
  pos_probs<-exp(forward(positive_hmm,test_observation))
  neg_probs<-exp(forward(negative_hmm,test_observation))
  
  pos_seq_prob<-sum(pos_probs[,final_index])
  neg_seq_prob<-sum(neg_probs[,final_index])
  
  if (pos_seq_prob>neg_seq_prob) incorrect<-incorrect+1
}

#最后计算交叉检验后的准确度
(cross_validation_accuracy <- 1 - (incorrect/nrow(promoters)))

例子2: 分析英语单词的字母顺序模式（寻找隐藏序列的patten）。
每个英语单词，其字母的组成顺序都是有一定规律的。我们很少看到在一个单词中，有两个元音字母连在一起。我们是否可以找到一个模型，来体现出来这种规律呢？
我们还是使用我们在上一讲中用到的文本数据。

#分析英语单词的字母顺序模式（寻找隐藏序列的patten）
library(ggplot2)
library("tm")
#读取文件
nb_pos <- VCorpus(DirSource(path_to_pos_folder), readerControl = list(reader = reader(DirSource(path_to_pos_folder)), language = "en"))
nb_neg <- VCorpus(DirSource(path_to_neg_folder), readerControl = list(reader = reader(DirSource(path_to_neg_folder)), language = "en"))
nb_all <- c(nb_pos,nb_neg) #合并两个数据
nb_all <- tm_map(nb_all, content_transformer(tolower)) #全部转换为小写

#因为nb_all的格式不好处理，我们将nb_all中的文本信息提取出来，存入到texts变量中
texts <- sapply(1 : length(nb_all), function(x) nb_all[[x]])

#文本处理
texts<-sapply(texts, function(x) gsub("\\s","W", x)) #把空格用W代替
texts<-sapply(texts, function(x) gsub("[0-9]","N", x)) #把数字用N代替
texts<-sapply(texts, function(x) gsub("[[:punct:]]","P", x)) #把标点符号用P代替
texts<-sapply(texts, function(x) gsub("[^a-zWNP]","O", x)) #其他符号用O代替
texts
#因为texts比较大，我们只取前40条数据，并将其每个字母打散，存入到列表中
big_text_splits<-lapply(texts[1:40], function(x) strsplit(x, ""))
big_text_splits<-unlist(big_text_splits, use.names = F) #将列表转换为向量

#由于我们不知道单词构造到隐藏状态，因此这里强行安了3个状态
states <- c("s1", "s2", "s3")
numstates <- length(states)
symbols <- c(letters,"W","N","P","O") #观测变量就是26个英文字母加上我们前面定义的几个特殊字母
numsymbols <- length(symbols)

#初始化一个初始状态概率矩阵pi。因为我们的隐藏状态也是随意给的，因此这里的初始状态概率矩阵也可以随机给一个
set.seed(124124)
startingProbabilities <- matrix(runif(numstates), 1, numstates)
startingProbabilities<-sweep(startingProbabilities, 1, rowSums(startingProbabilities), FUN="/")

#初始化一个转移状态矩阵
set.seed(454235)
transitionProbabilities<-matrix(runif(numstates*numstates),numstates,numstates)
transitionProbabilities<-sweep(transitionProbabilities, 1, rowSums(transitionProbabilities), FUN="/")

#初始化一个观测状态转移概率矩阵（发射矩阵）
set.seed(923501)
emissionProbabilities<-matrix(runif(numstates*numsymbols),numstates,numsymbols)
emissionProbabilities<-sweep(emissionProbabilities, 1, rowSums(emissionProbabilities), FUN="/")

#初始化隐马尔可夫模型
library("HMM")
hmm <- initHMM(states, symbols, startProbs = startingProbabilities, transProbs = transitionProbabilities, emissionProbs = emissionProbabilities)
#使用Baum-Welch 算法来学习参数。
hmm_trained <- baumWelch(hmm, big_text_splits)

hmm_trained$hmm

#根据学习后的参数，我们来画出转移概率矩阵。转移概率矩阵体现了在一个单词中，每个字母出现的概率
#画出从状态S1到各字母之间的概率
p1 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[1,]), aes(x = names(hmm_trained$hmm$emissionProbs[1,]), y = hmm_trained$hmm$emissionProbs[1,]))
p1 <- p1 + geom_bar(stat="identity")
p1 <- p1 + ggtitle("Symbol Emission Probabilities for State 1")
p1 <- p1 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p1 <- p1 + xlab("State")  
p1 <- p1 + ylab("Emission Probability") 
p1

#画出从状态S2到各字母之间的概率 
p2 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[2,]), aes(x = names(hmm_trained$hmm$emissionProbs[2,]), y = hmm_trained$hmm$emissionProbs[2,]))
p2 <- p2 + geom_bar(stat="identity")
p2 <- p2 + ggtitle("Symbol Emission Probabilities for State 2")
p2 <- p2 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p2 <- p2 + xlab("State")  
p2 <- p2 + ylab("Emission Probability") 
p2

#画出从状态S3到各字母之间的概率
p3 <- ggplot(data=data.frame(hmm_trained$hmm$emissionProbs[3,]), aes(x = names(hmm_trained$hmm$emissionProbs[3,]), y = hmm_trained$hmm$emissionProbs[3,]))
p3 <- p3 + geom_bar(stat="identity")
p3 <- p3 + ggtitle("Symbol Emission Probabilities for State 3")
p3 <- p3 + theme(plot.title = element_text(lineheight=.8, face="bold", vjust=2))
p3 <- p3 + xlab("State")  
p3 <- p3 + ylab("Emission Probability") 
p3

#获得隐含状态转移概率矩阵
(trained_transition_probabilities<-hmm_trained$hmm$transProbs)

#模拟给定隐马尔可夫模型的状态和观测路径。即根据模型构造一个句子（单词）
set.seed(9898)
simHMM(hmm_trained$hmm, 30)

完结～

【参考文章】
机器学习——HMM(隐马尔可夫模型的基本概念)
HMM学习最佳范例一：介绍
 隐马尔可夫模型(HMM) - 1 - 基本概念
 通俗理解，隐马尔可夫模型的前前后后
 HMM-前向后向算法理解与实现
 【NLP】隐马尔可夫模型三个基本问题相关算法实现

机器学习（十）——概率图模型之隐马尔可夫模型

1. 概论

2. 评估模型

3. 解码问题

4. 学习问题

5.R代码的实现

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读