请各位读者深吸一口气……呼……
开始……


（一）

阿黄是大家敬爱的警官，他性格开朗，身体强壮，是大家心目中健康的典范。

但是，近一个月来阿黄的身体状况出现异常：情绪失控的状况时有发生。有时候忍不住放声大笑，有时候有时候愁眉不展，有时候老泪纵横，有时候勃然大怒……

如此变化无常的情绪失控是由什么引起的呢？据警队同事勇男描述，由于复习考试寝室不熄灯与多媒体作业的困扰，阿黄近日出现了失眠等症状；与此同时，阿黄近日登陆一个叫做“xiaonei网”的网站十分频繁。经医生进一步诊断，由于其他人也遇到同样的考试压力、作息不规律的情况而并未出现情绪失控；并且，其它登陆XIAONEI网的众多同学表现正常，因此可基本排除它们是情绪失控的原因。黄SIR的病情一度陷入僵局……
最近，阿黄的病情有了新的眉目：据一位对手相学与占卜术十分精通的小巫婆透露，阿黄曾经私下请她对自己的病情进行诊断。经过观察与分析终于有了重大发现：原来阿黄的病情正在被潜伏在他体内的三种侍神控制！他们是：修罗王、阿修罗、罗刹神。

据悉，这三种侍神是情绪积聚激化而形成的自然神灵，他们相克相生，是游离于个体意识之外的精神产物，可以对人的情绪起到支配作用。每一天，都会有一位侍神主宰阿黄的情绪。并且，不同的侍神会导致不同的情绪突然表现。然而，当前的科技水平无法帮助我们诊断，当前哪位侍神是主宰侍神；更糟的是，不同的侍神(3个)与不同的情绪（4种）并不存在显而易见的一一对应关系。
所以，乍看上去，阿黄的病情再次陷入僵局……

我们怎样才能把握阿黄情绪变化的规律？
我们怎样才能通过阿黄的情绪变化，推测他体内侍神的变化规律？

关键词：两类状态：
情绪状态（观察状态）：放声大笑，愁眉不展，老泪纵横，勃然大怒
侍神状态（隐状态）：修罗王，阿修罗，罗刹神

（二）

阿黄的病情引来了很多好心人的关心。这与阿黄真诚善良的品格不无关系。
关于侍神的特点，占卜师和很多好心人找来了许多珍贵资料。其中很多人经过一段时间的观察与记录后，在貌似毫无规律的数据背后，发现了侍神与情绪之间的内在规律！！他们在多次观测后，建立在大量数据基础上，表现出宏观的内在联系！
由于这些好心人大部分是TONGJI大学的人，所以，这种规律被称作统计规律。这些人被称为统计学家（orz太土了）…………

具体的规律被概括为：
1. 每天，哪位侍神主宰与前一天侍神是谁有很大关系！即：前一个侍神会影响下一个侍神出现的概率多少。
三个侍神，两两之间的转化的几率的大小，我们总结在一个对应表中： 

便于某些自称为数学家的人计算，我们习惯于写成矩阵形式：


2.每位侍神主宰时，所对应的情绪出现也有一定的规律：某种侍神出现时，情绪的出现有一定的规律性。比如，如果今天的侍神是修罗王，那么阿黄放声大笑的概率是不同侍神对应的不同情绪出现的几率，我们也归纳在一张表中：


便于计算，也可以写成矩阵的格式。


3．由于每天的侍神状态，情绪状态都取决于前一天，所以，只有知道最开始阿黄的情绪状态，即初始状态，才有可能来推导随后的日子里，阿黄体内的侍神和情绪的变化情况。

我们假定，阿黄体内的侍神首先出现的概率满足下面的表格
修罗王   阿修罗  罗刹神
〔0.63    0.17    0.20〕

至此，我们已经掌握了研究阿黄情绪变化规律的所有信息，它包括：
两类状态：侍神状态，情绪状态
三种关系：侍神的转换关系，侍神与情绪的关系，侍神的初始状态

便于计算，我们用数学语言来定义，便于今后计算。

                   修罗王 阿修罗 罗刹神
初始状态矩阵：Π＝〔0.64    0.17    0.20〕

   状态转移矩阵：
A =  

两态混合矩阵（描述侍神状态和情绪状态对应关系）


（三）


知道了关于阿黄的这些信息后，我们能做些什么呢？

一、估计“放声大笑”－“老泪纵横”－“勃然大怒”出现的概率：

“放声大笑”－“老泪纵横”－“勃然大怒”是一种非常危险的组合！如果某三天内阿黄出现了“放声大笑”－“老泪纵横”－“勃然大怒”的情绪变化，会引起严重精神损伤！！那么，这种情绪组合出现的概率是多少呢？

看起来比较麻烦，因为“放声大笑”可以对应三种侍神；其余两个也可以。三天，三种侍神，所有可能为3*3*3＝27种。

根据全概率公式：所以：
P（笑－泪－怒）＝P（笑－泪－怒 | 修罗王－修罗王－修罗王）＋P（笑－泪－怒 | 修罗王－修罗王－阿修罗）＋P（笑－泪－怒 | 修罗王－修罗王－罗刹神）＋………………＋P（笑－泪－怒 | 罗刹神－罗刹神－罗刹神） 
一共27项相加，就可以把所有的可能计算出来！
可以想像，这样的计算是灾难性的。

当然，在计算机计算时，可以用递归法简化计算，降低复杂度。
①我们来把每一天，阿黄的情绪状态串起来，多天就形成一个状态序列。其中第t天的状态就是Ykt

②在计算序列中某一中间状态的概率时，用所有可能到达该状态的路径之和表示。


比如在t=2时间，状态为阿修罗的概率可以用下面的路径计算：



最后的观察状态的部分概率表示，这些状态所经过的所有可能路径的概率。比如：





③用αt ( j ) 表示在时间t时 状态j的部分概率。计算方法如下：


αt ( j )= P( 观察情绪 | 侍神是j ) * P(在t时间所有到j的途径)


其中两项相乘中的第一项我们由两状态混合矩阵就可以得到：



而后一项，我们需要前一项的结果来确定。这也表现了HMM每个环节的状态是基于前一环节状态的特点。



④如果说，每一环节都依赖前一环节，那么最初的环节，也就是最初的状态怎么来计算呢？
很简单。初始状态Π就派上它的用场了。
比如：第一天阿黄放声大笑，那么：
a1（修罗王）＝0.63×0.6＝0.378
a1（阿修罗）＝0.17×0.25＝0.0425
a1（罗刹神）＝0.20×0.05＝0.01
归纳成数学式：
⑤我们说过，后一天侍神的情况是由前一天来决定的。现在，如果知道了at（j），那么，后面一天的情况at＋1（j）也应该知道了吧。怎么来算at＋1（j）呢？
太枯燥了～还是继续上面的例子吧
比如，t＝2时，情绪为老泪纵横，侍神为阿修罗，那么：
a2（阿修罗）＝P（阿修罗时老泪纵横的概率）×P(所有到达阿修罗的概率)
其中：
P(所有到达阿修罗的概率)
＝a1（修罗王）×P(修罗王→阿修罗)＋a1（阿修罗）×P（阿修罗→阿修罗）＋a1（罗刹神）×P(罗刹神→阿修罗)
＝0.378×0.25＋0.0425×0.125＋0.01×0.675

这样，知道a1，知道at和at＋1的关系，迭代啊迭代，迭迭代代无穷匮矣……最终会找到你需要的状态的概率的。
（而且，在编成时运用迭代可降低算法复杂度，不太懂，不多扯……）



最后，解答“放声大笑－老泪纵横－勃然大怒”的概率
每天的三个概率对应三位侍神
第一天：放声大笑
(0.63 × 0.6) = 0.37800002
(0.17 ×0.25) = 0.0425
(0.2 ×0.05) = 0.010000001

第二天：老泪纵横
(((0.37800002×0.5) + (0.0425*0.375) + (0.010000001*0.125)) * 0.15) = 0.03092813
(((0.37800002*0.25) + (0.0425*0.125) + (0.010000001*0.675)) * 0.25) = 0.026640628
(((0.37800002*0.25) + (0.0425*0.375) + (0.010000001*0.375)) * 0.35) = 0.039965626
第三天：勃然大怒
(((0.03092813*0.5) + (0.026640628*0.375) + (0.039965626*0.125)) * 0.05) = 0.0015225002
(((0.03092813*0.25) + (0.026640628*0.125) + (0.039965626*0.675)) * 0.25) = 0.009509727
(((0.03092813*0.25) + (0.026640628*0.375) + (0.039965626*0.375)) * 0.5) = 0.01635469
所以，最终所有可能加起来，“放声大笑－老泪纵横－勃然大怒”的概率为
0.0015225002＋0.009509727＋0.01635469= 0.027386917
黄SIR暂时可以放心。

（四）

还有什么应用呢？


二、由观测状态推测最大可能性的隐状态，


            即：由阿黄情绪变化推测体内侍神的变化


1.       还是用穷举法吧




总之，变化是有限地，概率是知道地，每步是可算地，大小是可比地……


算一算，比一比，总能找到最可能的一条路。


只是随着天数增长、侍神数增长等，计算的复杂度会呈指数增加，这一点很不利。



又枯燥了～再举个例子：


比如，某三天，阿黄十分不幸地出现了“放声大笑－老泪纵横－勃然大怒”的情绪变化。


我们十分想知道：是怎样的侍神组合，最可能导致这种情绪？


最佳组合要取：


MAX{P（笑－泪－怒 | 修罗王－修罗王－修罗王）,P（笑－泪－怒 | 修罗王－修罗王－阿修罗）,P（笑－泪－怒 | 修罗王－修罗王－罗刹神）,………………,P（笑－泪－怒 | 罗刹神－罗刹神－罗刹神）}


27个里面比出一个概率最大的来，搞定～



2.       程序计算，我们用递归方法降低计算复杂度：


我们知道，对应于阿黄的每一个情绪状态，侍神的变化只有一个最佳路径，也许是这样：




每一条 部分最优路径都对应一个关联概率--部分概率 （相当于前面的那个中间概率）。与前面不同 是最有可能到达该状态的一条路径的概率。


① 我们定义： δ(i,t)是所有序列中在t时刻以状态i终止的最大概率。当然它所对应那条路径就是部分最优路径。  δ(i,t)对于每个i,t都是存在的。这样我们就可以顺藤摸瓜找下去，在序列的最后一个状态找到整个序列的最优路径。


②那么，最初的状态（t＝1）的最优路径是什么？我们还是要依赖初始状态矩阵Π


这和上次的算法是一样的。


δ(修罗王,1)＝0.378


δ(阿修罗,1)＝0.0425


δ(罗刹神,1)＝0.01



③那么，对于时间为t时刻的中间状态，如何来找它的部分概率（即最佳路径）呢？


再举个具体例子，t时刻，到达X状态的路径可能有ABC三条。



由此图可以看出，到达X的最优路径是下面三条中的一条：



(状态序列), . . ., A, X                              


(状态序列), . . ., B, X


(状态序列), . . ., C, X


我们就要比较：


P（到达A的最佳路径）×P（A到达X的概率）


P（到达B的最佳路径）×P（B到达X的概率）


P（到达C的最佳路径）×P（C到达X的概率）


再乘以P（X（某种侍神）对应的观测状态（某种情绪））就可以算得状态概率


④那么，在t时刻对应某种观察状态的概率记为δt（i），那么


第一天，由于没有先导，只能直接利用两状态转换矩阵计算：




第t天：


δt（i）＝MAX｛δt－1（j）×P（j状态→i状态）×i状态下对应观察状态概率｝


数学公式为：




这样，又可以迭迭代代无穷匮矣了～




最后，我们再来计算黄SIR“放声大笑－老泪纵横－勃然大怒”的最可能侍神组合：


第一天：放声大笑


修罗王(0.63 * 0.6) = 0.37800002


阿修罗(0.17 * 0.25) = 0.0425


罗刹神(0.2 * 0.05) = 0.010000001


第二天：老泪纵横


修罗王max ((0.37800002*0.5), (0.0425*0.375), (0.010000001*0.125)) * 0.15 = 0.028350003


阿修罗max ((0.37800002*0.25), (0.0425*0.125), (0.010000001*0.675)) * 0.25 = 0.023625001


罗刹神max ((0.37800002*0.25), (0.0425*0.375), (0.010000001*0.375)) * 0.35 = 0.033075


第三天：勃然大怒


修罗王max ((0.028350003*0.5), (0.023625001*0.375), (0.033075*0.125)) * 0.05 = 0.000708750


阿修罗max ((0.028350003*0.25), (0.023625001*0.125), (0.033075*0.675)) * 0.25 = 0.00558140


罗刹神max ((0.028350003*0.25), (0.023625001*0.375), (0.033075*0.375)) * 0.5 = 0.006201562



可见，第一天，修罗王主宰阿黄最为可能；


第二天，由修罗王变为罗刹神，造成阿黄老泪纵横的可能最大；


而第三天，继续由罗刹神主宰阿黄，造成勃然大怒的可能最大


所以，对应“放声大笑－老泪纵横－勃然大怒”最可能的侍神组合为：修罗王－罗刹神－罗刹神




尾声


一、了解了这个故事，我们就粗浅地了解了隐马尔科夫模型的构成。它包括两类状态，三种关系：Π（初始状态）A（状态转移矩阵）B（两状态混合矩阵）



二、HMM模型的初步了解就到这里。其主要功能有三：


1.根据已知的HMM找出一个观察序列的概率。（计算某种情绪组合出现概率）


2.根据观察序列找到最有可能出现的隐状态序列 （由情绪组合推导侍神组合）


3.从观察序列中得出HMM （这是最难的HMM应用。也就是根据观察序列和其代表的隐状态，生成一个三元组HMM ( Π,A,B)。使这个三元组能够最好的描述我们所见的一个现象规律。限于水平，不讨论）



三、上面应用中，1设计到FORWARD算法，2设计Viterbi算法，可以查阅资料，都可以在程序中实现。

四、生物信息中涉及到HMM的应用有很多。在蛋白质DOMAIN描述，分子进化树的构建中都会应用这个模型。其实大同小异：我们能看得见的就是观察状态：序列（情绪），而隐藏在当前序列背后的各种隐状态：比如基因……（修罗王，阿修罗，罗刹神……）。基因发生了怎样的进化过程，我们不能直接观测，但可以用应用3――从观察序列中得到HMM，用已知序列进行训练，生成包含实际序列信息的HMM模型，就可以通过已有的序列来进行基因的人工注释（大概是这个意思，FT）



五、感谢来自这个网站：


www.comp.leeds.ac.uk/roger/HiddenMarkovModels/html_dev/main.html


提供的大量介绍，数据，图表，小程序，很有帮助


特别感谢丁香园里的崔晓源进行的详尽翻译和解释，这篇文章的思路归功于他的翻译作品。



六、感谢主人公姓名雷同者的大力支持，体谅和鼓励，以及所有带来灵感的同学～你们让生活更精彩！