非常崩溃！多序列前-后向算法的HMM训练问题？

你给的公式不是多序列比对的，是一个序列对的比对，只是数据有多个序列对，这种不叫做多序列比对。 Scaling是为了防止出现INF的情况，可以先不用scaling的方法试一试。不用scaling，一般也不怎么会出现underflow的情况，但是要用对数哦。 最后P_k为什么要用log后，再用指数呢？好像没必要。 序列比对，推荐一本书，Biological sequence analysis: probabilistic models of proteins and nucleic acids。 多序列比对一般是先两两比对，自底向上融合结果，最后做成multiple sequence alignment的。

HMM的EM学习算法的思路：假设知道了序列的比对，直接统计比对的各种state和emission即可。HMM模型给每一种比对一个概率，使用概率进行统计。 另外，建议要明确模型的状态和状态的emission。state和emission不同，hmm的实现需要具体设计。

EM训练的时候，Q函数的增量应该是满足李普希茨条件的，如果训练时，Q函数增量的序列不满足李普希茨条件，多半是模型有问题或者代码写得有问题，要耐心查，不能急。

首先严重感谢楼上大牛的回答！ 按你的意思，我理解的多序列出现了偏差。我想表达的是同一个HMM模型，独立的跑出K个序列样本： O11 O12 .... O1T O21 O22 .... O2T .... OK1 OK2 .... OKT 而要求是，通过上述的样本进行学习，并归纳出状态个数。这里的状态应该是指隐状态个数吧？按我的理解，如果在知道状态个数的情况下应该是通过文档里最后给的式子来重估HMM模型吧？ 楼上说的“序列的比对”是用于哪方面的？是HMM模型初始化吗？ 还有，文档里有一词叫 "multi stream"又是指什么？

我主要做基因组的序列比对算法。 在DNA序列分析里，序列是一串A、C、G、T字符。比如说有两条序列，序列a=AACGTT，序列b=ACGG，序列比对是找出a和b的字符一一配对的最佳方式，允许在序列中插入空字符‘-’。例如，a和b的一种比对方式， a: AACG-TT b: A-CGG-- 序列比对分为序列对比对（Pair-wise alignment）和多序列比对（Multiple alignment）。上面的例子是序列对比对。多序列比对的例子见链接http://en.wikipedia.org/wiki/Multiple_sequence_alignment 序列比对算法的两个核心部分，其一是打分，其二是动态规划算法。打分，英文叫做scoring function，因为两条序列的可能的比对很多，打分策略给每种比对一个分数，分数最高的作为保留为最后的结果。不同的打分策略，挑出来的比对不同。动态规划算法，经典的序列比对算法Needleman-Wunsch和Smith-Waterman算法都是基于动态规划，我不清楚有没有序列比对算法是不用动态规划来做的。 序列对比对的具体实现是在一个二维的动态规划表上寻找最优比对。多序列的比对则是在一个N维的动态规划表上进行寻优，复杂度是O（N^d），N是序列的长度，d是序列的数目。这基本上是无法做的事情，所以会用我之前说的方法，先做序列对比对，把比对后的结果再进行配对进行下一轮的序列对比对，这样反复，最后得到多序列比对的结果。有点Hierarchical clustering的感觉。 你说的例子，有点Profile-HMM的味道，是一种序列对比对的变异体。当然，这只是我的感觉，因为不清楚你的具体背景，所以有可能感觉错了。我就照着我的感觉走了。做HMM最重要的一步是要明白自己的模型是什么，模型包含多少个状态（就是你说的隐状态），每个状态的行为是什么。这些明确了，HMM也就做好了49%，剩下就是coding，占1%吧，debug占50%，debug实在太痛苦了。模型的逻辑对了，写出来的代码一般不会有打错。所以理清自己的逻辑非常重要。 HMM的viterbi算法，forward算法，backward算法，forward-backward算法都是动态规划，不会很难。

再来说说，HMM的训练。 初始点的选择，可以使用动态规划计算出一种比对作为初始。 HMM的EM训练算法，可能之前我描述的不是很清楚，就以你给的例子为例吧，比如给定一个序列a=a1a2...an，然后通过HMM产生了K个样本 O11 O12 .... O1T O21 O22 .... O2T .... OK1 OK2 .... OKT 现在有一个非常聪明的算法能够告诉你每个样本是怎样通过HMM产生的，也就是相当于知道“真实”的比对 1）O11 O12 .... O1T a1 a2 ...........an 2) O21 O22 .... O2T - a1...... an .... K) OK1 OK2 .... OKT a1 - .......-.....an 知道了这些“真实”的比对后，HMM参数估计就是一个简单的计数，然后归一化的操作。可是，问题是，“真实”的比对不知道，所以会用EM算法来学习参数。因为EM能够给出一个序列对的每一种可能的比对的概率，因此，就可以去数每种状态转移的期望次数，每个状态的每种emission的期望次数，然后归一化。不断的迭代直到收敛。 我刚开始做的时候，花了两个月写了最基本的动态规划，经常出错，不停的在反复的琢磨逻辑是不是正确。熟悉了动态规划后，写HMM就很快，但是和你现在的情况很像，Q函数的delta忽大忽小。慢慢的清楚了，基本上是自己的逻辑没有理清楚，有时候是因为代码上一些不太容易看出的细节。 我的建议，明确自己的问题，想清楚模型，不要急着写代码，把模型的每个细节都考虑清楚，我自己用了很笨的办法，现在纸上写一遍代码，自己觉得没错了再code。 如果，你是做工程赶时间的话，我建议直接找找现成的适合你的问题的代码，应该网上都能找到。

如果是在做research，不在乎时间就是为了弄清楚HMM，建议在中间加入print out，把每一步的transition matrix和emission matrix输出，check每一步结果是不是合理。

对不起，看到序列，我自然反应成了序列比对了，呵呵。我看了一下Rabiner的tutorial，可以理解你说的多序列是什么了。Rabiner在描述算法的时候，用了只有一个样本的例子。在后面的推广，引入了Multiple observed sequence的概念，你理解的多序列没错，就是一个HMM产生了一堆数据。你给的文件里面红色框框的公式也没错，公式的意思就是我上面说的意思，只是我比较直白。

其实我在做HMM完全只是帮别人解决麻烦，自己也不是这方面专业的，因此，只想着快速得到想要的结果。 网上确实有一些代码，最近理了一下思路，再结合别人的代码，倒也有了一些结果，只是HMM估计的的值不太准确，比如2个隐状态，2个观察值下， T=500,K=20的观察值用于测试，重估的结果也不算太准，A矩阵差别较大， 比如（0.5， 0.5）的概率估出来可能是（0.2，0.8），不知道结果算不算合理？ 况且我看题目的要求还挺奇怪的，说是要从多序列里归纳有状态个数，但我看过的书都是认为状态个数是预先知道的; 还有需要用学习好的HMM模型来预测某序列的下一个状态是什么？这点我就更难理解了：重估出来的是一个概率矩阵，这怎么可能去“预测”下一个状态？应该只能说下一状态以多少概率发生吧？ 楼上应该是北邮生物电子专业的？