设为首页收藏本站
开启辅助访问 切换到宽版

科学网

 找回密码
  注册
科学网 › 标签 › HMM

tag 标签: HMM

相关帖子

 版块 作者 回复/查看 最后发表

没有相关内容

相关日志

隐马尔科夫模型简介(三)
luria 2019-7-15 20:33
本篇继续讨论隐马尔科夫模型的第二个基本问题——解码问题 在讨论之前，还是欢迎我们的主角——隐马尔科夫模型 (HMM) 出场！ Problem 2. 解码问题 ( Decoding Problem) 已知 HMM 模型λ =(A, B, π ) ，以及观测序列 O ，求最可能的状态序列。 之前的博文 http://blog.sciencenet.cn/blog-2970729-1188964.html 中讨论的问题就属于这一类问题。 其实这个问题和之前讨论的评估问题 (第一类问题) 相比，已知的条件是一致的，只是所求不同。因此理论上也可以用枚举法，事实上，在第一篇 隐马尔科夫模型 简介博文中我们也这样做了。同样因为计算量比较大，所以枚举法并不可取。在此介绍一种目前主流的解决方案，叫做 维特比算法 (Viterbi Algorithm) ，算法演示如下： Viterbi Algorithm.gif Viterbi 算法的思想是步步为营，每一步都是当前最佳的路。如果用气温和年轮的例子来说： 1. 首先写出第一年的情况 因为第一年观测到的是小年轮，于是 δ 1 (H) 表示第一年观测到小年轮的同时气温为高温的概率： δ 1 (H) = 0.6 × 0.1=0.06 注：其中 0.6 表示初始为 H 的概率， 0.1 表示观测到小年轮而且为高温的概率。 同时，第一年观测到小年轮而且气温为低温的概率为： δ 1 (C) = 0.4 × 0.7=0.28 δ 1 (C) δ 1 (H) ，因此 第一年为低温的可能性更大 ，所以下一步计算时不考虑第一年为高温的情况 2. 计算第二年的情况 第二年观测到的是中年轮 (M) ，于是 δ 2 (H) = δ 1 (C) × 0.4 × 0.4 = 0.0448 注：第一个 0.4 表示由低温转高温的概率，第二个 0.4 为观测为中年轮，且状态为高温的概率 同理，δ 2 (C) = δ 1 (C) × 0.6 × 0.2 = 0.0336 δ 2 (H) δ 2 (C) ，因此 第二年为高温的可能性更大 ，所以下一步计算时不考虑第二年为低温的情况。 3. 计算第三年的情况 第三年观测为小年轮，计算方法同第二年，直接写出结果如下： δ 3 (H) = δ 2 (H) × 0.7 × 0.1 = 0.003136 δ 3 (C) = δ 2 (H) × 0.3 × 0.7 = 0.009408 δ 3 (C) δ 3 (H) ，因此 第三年是低温的可能更高 ，所以下一步计算时不考虑第三年为高温的情况。 4. 计算第四年的情况 第四年观测为大年轮，计算方法同第二年，直接写出结果如下： δ 4 (H) = δ 3 (C) × 0.4 × 0.5 = 0.0018816 δ 4 (C) = δ 3 (C) × 0.6 × 0.1 = 0.00056448 δ 3 (H) δ 3 (C) ，因此 第四年是高温的可能更高 。综上，这四年最可能的状态是 C-H-C-H ，这一结果和第一篇隐马尔科夫简介博文( http://blog.sciencenet.cn/blog-2970729-1188964.html )中 HMM 算法得到的结果一致。但是通过计算过程来看，更加高效。 如果用数学公式表示，如下： 参考材料： https://sens.tistory.com/m/320?category=501289 Guy Leonard Kouemou. History and Theoretical Basics of Hidden Markov Models Mark Stamp. A Revealing Introduction to Hidden Markov Models 李航 . 统计学习方法
个人分类: Algorithm|2223 次阅读|0 个评论
隐马尔科夫模型简介(二)
luria 2019-7-15 20:21
本文接上篇继续讨论隐马尔科夫模型，上篇详见： http://blog.sciencenet.cn/blog-2970729-1188964.html 在讨论之前，欢迎我们的主角——隐马尔科夫模型 (HMM) 出场！ 本篇主要探讨隐马尔科夫模型三个基本问题中的第一个问题，后续还会接着讨论另外的两个，敬请关注 3. 隐马尔科夫模型的三个基本问题 Problem 1. 评估问题 ( Evaluation Problem) 给定模型λ = (A, B, π ) ，以及观测链 O ，求解此模型出现观测链 O 的概率 P(O| λ )? （注：用上篇中年轮与气温的例子来说，已知状态转移矩阵 A ，观测矩阵 B ，以及初始状态分布π，求特定的观测链，例如 S-M-S-L 这种年轮排列，出现的概率） 方案 1. 枚举法 为了解决这个问题，我们最开始可能想到的是简单粗暴的枚举法 ( 概率法 ) ：列举出所有可能的状态链 ( 例如 H-H-H-H 、 H-H-H-C 、 H-H-C-H … 等 ) ，对应出现 S-M-S-L 这种年轮排列的概率，再把所有的概率相加，自然就可以得到在 P(O| λ ) 。 如果按照数学公式来计算，思路如下： 先定义状态链 X = (X 0 , X 1 , X 2 , …, X T-1 ) 。 (1) 只看 HMM 图中虚线上半部分，相当于给定模型λ，求状态链 X 的概率，很显然 (2) 只看 HMM 图中虚线下半部分，相当于在给定了状态链 X 和模型λ，求观测链 O 的概率。很显然 (3) 按联合概率的定义有 当然，这个公式也可以推导出来，具体如下： (4) 再按概率公式展开 P(O| λ ) ，即 也即， 这种算法虽然直观，但是计算量非常大，是 O(TN T ) 阶 。在概念上可行，但是计算不上可行。 方案 2. 前向算法 (Forward Algorithm) 既然叫做隐马尔科夫链，作为一条链我们就可以用链的方式来算。链的特征是当前的状态仅与上一个状态有关。 Q1 ：还是用年轮与气温的例子来说，如何计算第四年观测到大年轮 (L) 的概率？ A1 ：它是以下两部分之和：其一是第四年气温为高温 (H) ，观测为大年轮 (L) 的概率，其二是第四年气温为低温 (C) ，观测为大年轮 (L) 的概率。 Q2 ：这时有人要问了，第四年气温为低温 (C) 的概率是多少？ A2 ：它是以下两部分之和：其一是第三年气温为 H 的概率乘以气温由高温 (H) 转低温 (C) 的概率，其二是第三年气温为低温 (C) 的概率乘以气温由低温 (C) 转低温 (C) 的概率。 Q3 ：紧接着问，第三年气温为高温 (H) 的概率是多少？ A3 ：这是第三年观测为小年轮 (S) ，气温为 H 的概率！ Q4 ：那么第三年观测为小年轮 (S) 的概率又是多少？ A4 ：很好，回头看看第一问 Q1 ，这个问题又回来了！相当于整个计算递归 ( 轮回 ) 了一遍。 用图来说明如下： 如果用数学公式表示，如下： 先定义 前向概率 α t ( i ) ，它表示到时刻 t 时，观测序列为 o 1 , o 2 , o 3 , …, o t ， 且状态为 q i 的概率，记作： 思路采用的是递归 ( 从最后一个状态到第一个状态 ) ，但是计算的时候反过 ( 从最一个状态到最后一个状态 ) ： 1. 计算初始的前向概率α 0 (i) 2. 联系 t-1 时刻的前向概率 α t-1 (j) 与 t 时刻的前向概率α t (i) 注：它由两部分组成，红色框中的部分表明 t-1 时刻各个状态对应的 前向概率与其状态转移概率之积再求和，前面 t-1 的事就搞定了；蓝色框内的是当前时刻 (t 时刻 ) O t 这种观测值的第 i 个状态的概率 3. 当运行到最后一个前向概率时，所有的计算都完成了，即获取到 P(O| λ ) ，具体的式子如下： 直接概率法的复杂度是 O(TN T ) ，而前向算法的复杂度是 O(N 2 T) ，当 T 较大时可以大大节省计算开销。 参考材料： Guy Leonard Kouemou. History and Theoretical Basics of Hidden Markov Models Mark Stamp. A Revealing Introduction to Hidden Markov Models 李航 . 统计学习方法
个人分类: Algorithm|2614 次阅读|0 个评论
隐马尔科夫模型简介(一)
luria 2019-7-10 18:43
1. 一个简单的例子 假设我们想知道某个固定的地区一些年来的平均年平均气温。为了让这个问题变得有趣，我们假设我们关注的这些年份是在温度计发明之前的遥远过去。我们不能回到过去，所以需要寻找与温度相关的非直接证据。 为了简化问题，仅会考虑两种年平均温度， hot 和 cold 。假设现代的证据证明 hot year 之后又是一个 hot year 的概率是 0.7 ， cold year 之后又是一个 cold year 的概率是 0.6 ，而且这个概率在过去的岁月里都适用，那么目前的信息总结如下： （ 1 ） 其中 H 代表 hot ， C 代表 ''cold 。 另外，当前的研究证实树的年轮与温度相关。也为了简化起见，我们将树的年轮分为三种大小： small ， medium 和 large ，分别用 S, M 和 L 表示。最后基于存在的证据显示，年平均温度与树的年轮之间的关系如下： （ 2 ） 对于这个例子来说， 状态 (state) 是年平均气温， H 或 C 。从一种状态到另一种状态的转移过程是 马尔科夫过程 (Markov process) 。因为下一个状态仅依赖于当前状态，而且符合如矩阵 (1) 的固定概率。然而，真实的状态是隐藏的 (hidden) ，因为我们不能直接获取到过去那些年份的气温。 虽然我们不能观测到过去那些年的状态（气温），但是我们可以观测到树的年轮大小。通过矩阵 (2) ，树的年轮告诉我们关于气温的概率信息。因为状态是隐藏的，这种类型的系统我们称为 隐马尔科夫模型 (Hidden Markov Model ， HMM) 。我们的目标是有效地，且高效地利用观测到的数据了解马尔科夫过程的不同特征。 我们称矩阵 (1) 为 状态转移矩阵 (state transition matrix) ： （ 3 ） 称矩阵 (2) 为 观测矩阵 (state transition matrix) ： （ 4 ） 在这个例子中，还需要知道一个 初始状态分布 (initial state distribution) ，记为 π，而且假设 （ 5 ） 这里的初始状态分布指关注的这几年中，最开始的那一年，高温的概率为 0.6 ，低温的概率是 0.4 ( 注：这个对应关系在文章中虽未明说，但从后面的算式中可查看到) A 、 B 和π中每个元素都是概率值，矩阵中每行都是一个概率分布，每行之和为 1 。 现在我们关注的四年 ( 例如 2007-2010 年 ) ，我们观测到这四年树的年轮分别为 S, M, S 和 L ，且用 0 表示 S ， 1 表示 M ， 2 表示 L ，那么观测链如下： （ 6 ） 如下图： 通过观测到的年轮结果，我们想推测出最可能（ most likely ）的马尔科夫过程状态链 ( 注：也即这四年气温情况分别是怎样的 ) 。这个问题不像看起来那么可以一刀切（ clear-cut ），因为对 most likely 人们有不同的理解。一方面，我们可以将 most likely 大致地定义为，所有长度为 4 的可能状态链中概率最高的那种状态链。 ( 注：也即所有可能的 2^4 种状态 ，每种状态都有对应的概率，这其中概率最大的那种 ) 。 我们可以采用动态规划算法 (Dynamic programming, DP) 来求解出这个特殊解。 另一方面，我们可能会将 most likely 合理地定义为，使状态链中每个状态的可能最大化 ( 注：例如状态链中第一个状态最可能是低温，第二、三、四个状态最可能是高温、低温、高温，那么最后整条状态链最可能是低温 - 高温 - 低温 - 高温 ) 。 其实，动态规划算法得到的解与隐马尔科夫算法得到的解并一定完全一致。动态规划算法必须要有合理的状态转移，而这个在隐马尔科夫模型算法中不是必须的。 ( 注：老实说读到这里才开始真正佩服 HMM ！另外，有关动态规划算法，欢迎关注我的另一篇博文： http://blog.sciencenet.cn/blog-2970729-1112311.html ) 而且，即使所有的状态转移都是合理的， HMM 方法也不同于 DP 方法，那就是在处理 标记问题 时（ notation ）。 2. 标记问题 所有可能的观测结果通常用｛ 0,1,...,M-1 ｝表示，因为这样可以在不失普遍性的前提下， 使标记问题简单化 。观测到的序列依次用 O i 表示，其中 O i ∈ V for i = 0, 1,...,T-1 。 注：用上面的例子来说， T 为 4 ，表示观测的年份数，从 2007-2012 共 4 年； N 为 2 ，表示观测状态数，气温要么是高温，要么是低温； M 是 3 ，表示树的年轮大小数，年轮要么是大号 (L) ，要么是中号 (M) ，要么是小号 (S) ； Q 是高温和低温，即 {H ， C} ，它是一个集合！是所有可能的马尔科夫状态集合； V 是大号 (L) 、中号 (M) 和小号 (S) ，对应为数字即 {0, 1, 2} ，是一个集合！是所有可能的观测结果集合； A 是状态转移矩阵，是一个矩阵！如果写成概率公式，如下： B 是观测矩阵，是一个矩阵！如果写成概率公式，如下： 注意：这两个公式后面会有用到！ 另外， A 和 B 矩阵是行随机的。它们和状态无关，是客观存在的规律。 π是初始状态分布； O 是观测的序列，例如 (S, M, S, L) ，或者换算为数字 (0, 1, 0, 2) 整个马尔科夫过程用图展示如下： 其中 X i 代表隐状态，其它的解释如上。隐藏于虚线上方的状态链，即马尔科夫过程，取决于当前状态和 A 矩阵。我们仅能观测到 Oi ，它通过 B 矩阵与隐状态关联起来。 HMM 可以通过 A, B 和 π确定，因此可以用以下公式表示： λ = (A, B, π ) 那么，一个长度为 4 的状态链可以写为： 例如：低温，高温，低温，高温。这个是我们最终想要知道的。 与其相对应的观测链如下： 例如：年轮 S, 年轮 M, 年轮 S ，年轮 L 。是我们观测到的 于是，状态 x 0 的概率为π x0 。在 x0 的条件下观测到 O 0 的概率是 b x 0 ( O 0 ) ， 从 x0 转移到 x 1 的概率为 a x0,x1 。依次类推，最终观测为 O ，推测为 X 这种链 ( 例如高温 - 高温 - 低温 - 低温 ) 的概率是： 在上面通过年轮推测气温的例子中，观测 O=(0,1,0,2) ，对应 ( 小年轮，中年轮，小年轮，大年轮 ) 。那么高温 - 高温 - 低温 - 低温这种可能的状态链的概率是： P(HHCC) = 0.6(0.1)(0.7)(0.4)(0.3)(0.7)(0.6)(0.1) = 0.000212 同理，我们可以算出每种可能情况 ( 总的是 2^4 种可能，有四年，每年两种可能的温度情况 ) ，如下表： 其中， normalized probability 是这 16 种 可能的状态链的概率归一化之后的结果。 如果是动态规划算法， 推测最可能状态链是 CCCH ，因为其 normalized probability 最大；如果是隐马尔科夫模型，则计算如下： 1. 先算出上表中所有第一个状态为 H 的可能性之和，即 HHHH ， HHHC ， HHCH ， HHCC ， HCHH ， HCHC ， HCCH ， HCCC 的 normalized probability 之和，为 0.188182 。第一状态为 C 的可能性之和为 1-0.188182=0.811818 2.再依次算出第二个状态，第三个状态和第四个状态的 H/C 概率，如下表： 因为链的各个 element 计算是独立的，这里只需要挑选出每个 element 概率最大的状态，构成概率最高的状态链，本例中 HMM 计算出的最可能状态链是 CHCH 。 注意：两种算法计算出的最可能状态链不同，但是两种状态链都是合理的。 参考文献： Mark Stamp. A Revealing Introduction to Hidden Markov Models
个人分类: Algorithm|7382 次阅读|0 个评论
单荧光分子光漂白数据分析
grating 2014-8-24 17:54
5120 次阅读|0 个评论
后及性、无前溯性与无后效性——马尔可夫过程“无后效性”回溯谈
hillside 2013-2-27 10:58
后及性、无前溯性是法律等方面涉及时间效力的专业用语，而“无后效性”是统计学马尔可夫过程（或称马尔可夫链）的用语。 我们来看“无后效性”比较通用的一个定义：如果一个过程的“将来”仅依赖“现在”而不依赖“过去”，则此过程具有马尔可夫性,或称此过程为马尔可夫过程。 我总觉得“无后效性”看上去有些别扭，不如“无前溯性”来得直白与爽快，因为“后效”往往正是我们需要的或者难免避免的“下一次”。文革期间，“以观后效”类用语曾盛极一时。 “无后效性”的抽象与假设使马尔可夫过程取得了很大的成功。不过，成也萧何，败也萧何。割断历史、无视历史的代价也是显著的。 隐马尔可夫过程对马尔可夫过程进行了补救。 附：下面转引网络文章《爱情的隐式马尔可夫模型(Love in the Hidden Markov Model)》片断： 爱情的隐式马尔可夫模型(Love in the Hidden Markov Model) 首先感谢原英文作者Tom Yeh的精彩描述，生动地讲述了HMM模型的原理,在此我斗胆用我自己的语言用中文修改描述一次。 男生和女生分别是来自不同星球的科学事实已经众所周知的了。男生们总是认为,女生们都是谜一样的生物，他们的情感状态浮动似乎是以秒单位在变化的，难以理解，更勿论预测了! 而女生们觉得男生都是没有感觉动物,完全不能理解什么叫感受-尽管已经告诉他们N次了!这种男女之间的根本差别，导致了他们之间的感情关系是受一种超级无敌复杂的系统所支配的。 不过,我们可以用一个叫隐式马尔可夫(Hidden Markov Model)的数学模型来分析这个系统。 小明,作为一名计算机科学家, 决定要系统地去分析他女朋友的情感不确定性, 挖掘出里面的规律!于是乎,小明仔细地记录了半年来他女朋友小丽每天的喜怒哀乐变化状态, 并作了一张图表来表示小丽的历史情感变化。小明想知道, 有了这些数据,他能否从中得出知道, 如果小丽某天的情感状态是高兴, 那么第二天她更多的是保持好心情呢，还是更多地变得悲伤了.如此等等... 数据胜于雄辩。小明从这半年的数据里面发现，当小丽高兴的时候，3/4的情况下第二天她仍然保持着好心情，只有1/4的情况小丽第二天心情会改变，比如变得气愤,悲伤等等。小明继续分析其他各种情感状态变化情况，比如从高兴到悲伤, 悲伤到气愤, 高兴到气愤等所有的可能组合。很快小明就得到所有的组合变化数据,从中得知对于任意小丽的某天情感状态下,下一个最有可能的情感状态。 …… 这个过程,同学们,就是有名的 马尔可夫过程 (Markov process) 不过需要注意的是, 马尔可夫过程有一些假设的前提. 在我们的例子里面, 预测下一天小丽的心情, 我们只依赖当天小丽的心情,而没有去考虑更先前她的心情。很明显这种假设下的模型是远不够精确的。很多时候,随着日子一天一天的过去,女生一般会变得越来越体谅。经常女生生气了几天后,气就会慢慢消了. 比方说如果小丽已经生气了3天了,那么她第二天变得高兴起来的可能性,在多数情况下,要比她只生气了一天而第二天变得高兴的可能性要高. 马尔可夫过程并没有考虑这个, 用行话讲, 就是马尔可夫模型忽略远距离历史效应 ( long range dependency)。 有些时候，我们无法直接观测一个事物的状态。比方说, 有些女生是很能隐瞒自己的情感而不流露出来的! 他们可能天天面带微笑但不代表他们就天天高兴。因此我们必须要有窍门, 去依赖某些我们能够直接观察到的东西。 话说回来，我们的主人公小明，自从被小丽发现他这种近乎变态的科学分析行为后，变得非常善于隐藏自己的心情，导致某天小明错误估计了小丽的心情!在误以为那天小丽会心情好的情况下，小明告诉小丽自己不小心摔坏了她心爱的iPod...，小明没想到其实那天小丽正因为前一天错过了商场名牌打折扣的活动而异常气愤... 一场血雨腥风过后,两个人最终分手了。 不过，小明凭着自身的英俊高大潇洒，很快又交上了另外一个女朋友 - 小玲。鉴于小明意识到，女生表面的情感流露非常不可靠，小明决定要另寻他径, 继续预测女朋友的心情! (作为一个科学家，小明的确有着不怕碰壁的精神!)小明每个月都帮小玲付信用卡的费用(真不明白，有这样的男朋友，小玲有什么理由不高兴啊!), 因此小明每天都可以通过Online banking知道小玲每天都买了什么东西。小明突然灵机一动: 没准我能通过观测她的购物规律，推导预测出小玲的心情! 听起来有点匪夷所思，不过这个过程，的的确确是可以使用叫作隐式马尔可夫的数学模型来表示并分析的。 由于我们需要预测的变量 - 心情状态是无法直接观测的，是隐藏 (Hidden) 起来的，因此这种模型才叫隐式马尔可夫模型。 隐马尔可夫模型在计算机语音识别等领域取得了惊人的成功。据称：“到目前为止，HMM(隐马尔可夫模型)一直被认为是实现快速精确的语音系统的最成功的方法。” 在地球科学领域， 隐马尔可夫模型在降水、气温等大气因素的预测上已经有了一些探索。在水文学方面的应用还很薄弱，相信研究会逐步增多。用地球科学特别是近年火爆的气候变化研究术语来看，马尔可夫过程的“明”与“隐”可以比作“代用指标”与“本征指标”(注：我一时不知对应“代用指标”的原始指标的通用说法是什么，但我觉得“本征指标”挺妙)。
个人分类: 数学与统计园地|16285 次阅读|0 个评论
基于HMM的语音合成技术中的参数生成算法
naxingyu 2013-1-14 20:26
一、定义与定理 在基于隐马模型的语音合成技术中，连续密度隐马尔科夫模型（CD-HMM）集用于将语音参数建模，每个HMM状态的输出状态用单高斯函数（Gaussian）或混合高斯函数（GMM）表示(Zen et al., 2009)，其参数生成算法的目标是在给定高斯分布序列的前提下，计算出具有最大似然函数的语音参数序列(Tokuda et al., 1995)。$$p(\mathbf{o}\mid\lambda)=\sum_{all\mathbf{q}}p(\mathbf{o}\mid\mathbf{q},\lambda)P(\mathbf{q}\mid\lambda)$$其中$\mathbf{o}:=\{\mathbf{o}_1,\mathbf{o}_2,\ldots,\mathbf{o}_T\}$代表语音参数矢量序列，$\mathbf{q}:=\{q_1,q_2,\ldots,q_T\}$表示高斯分布序列。本文中，大写粗体字母表示矩阵，如$\mathbf{W}$，小写粗体字母表示矢量，如$\mathbf{o}$，普通字母表示标量。小写字母$p$专用于表示连续变量的概率分布，大写字母$P$专用于表示离散变量的概率分布。除非特别说明，所有的矢量都是列矢量。符号$\top$用于表示矩阵或矢量的转置，如$\mathbf{W}^{\top}$。由均值矢量$\boldsymbol{\mu}$和方差矩阵$\boldsymbol{\Sigma}$表示的生成数据点$\mathbf{o}_t$的高斯密度函数写作$\mathcal{N}(\mathbf{o}_t\mid\boldsymbol{\mu},\boldsymbol{\Sigma})$。 在传统HMM中，状态（高斯分布）序列是由转移概率矩阵决定的，既 $$P(\mathbf{q}\mid\lambda)=\prod_{t=1}^TP(q_t\mid q_{t-1},\lambda),$$ 在基于HMM的语音合成中，状态序列是由显式时长模型输出的时长特征矢量决定的。由于这种设定改变了模型的严格马尔科夫性，我们将其成为隐半马尔科夫模型（HSMM）。因此，下面的推导专注于输出概率密度函数。 在典型的语音识别和语音合成系统中，声学参数按帧提取，第$t$帧的参数由矢量表示为 $$\mathbf{c}_t=\left ^{\top},$$ 第$t$帧的观测值（对于模型来说，即输出值）通常定义为由声学特征及其一阶和二阶动态特征共同组成的矢量 $$\mathbf{o}_t=\left ^{\top}.$$ 这些动态特征是以相邻帧静态特征的回归系数的形式计算得到的，即 \begin{aligned} \Delta\mathbf{c}_t = \sum_{\tau=-L_-^{(1)}}^{L_+^{(1)}}w^{(1)}(\tau)\mathbf{c}_{t+\tau},\\ \Delta^2\mathbf{c}_t = \sum_{\tau=-L_-^{(2)}}^{L_+^{(2)}}w^{(2)}(\tau)\mathbf{c}_{t+\tau}. \end{aligned} 因此，特征参数序列可以表示为 \begin{aligned} \mathbf{o} = \left ^{\top},\\ \mathbf{c} = \left ^{\top}. \end{aligned} 观测值$\mathbf{o}$和静态特征$\mathbf{c}$之间的关系为 $$\mathbf{o}=\mathbf{W}\mathbf{c}$$ 其中$\mathbf{W}$定义为窗系数矩阵，即 \begin{aligned} \mathbf{W} = \left ^{\top}\otimes\mathbf{I}_{M\times M},\\ \mathbf{W}_t = \left . \end{aligned} $\mathbf{w}_t^{(d)}$是用于计算第$t$帧第$d$阶动态特征的窗系数，只在第$t$位和相邻位置的元素有非零值，非零值的范围取决于窗宽度，通常为1，即前后一帧。 $\mathbf{I}_{M\times M}$表示$M\times M$单位矩阵，用于将相同的窗系数应用于所有$M$维参数。 因此，似然函数表示为 $$p(\mathbf{o}\mid\mathbf{q},\lambda)= \prod_{t=1}^{T}p(\mathbf{o}_t\mid q_t,\lambda)= \prod_{t=1}^{T}\mathcal{N}(\mathbf{o}_t\mid\boldsymbol{\mu}_{q_t},\boldsymbol{\Sigma}_{q_t})= \mathcal{N}(\mathbf{o}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}})$$ 其中 \begin{aligned} \boldsymbol{\mu}_{\mathbf{q}} = \left ^{\top},\\ \boldsymbol{\Sigma}_{\mathbf{q}} = diag\left ^{\top} \end{aligned} 根据以上内容，我们有如下定义 \begin{aligned} \text{static feature} ~~~~ \mathbf{c} ~~~ MT\times 1\\ \text{observation} ~~~~ \mathbf{o} ~~~ 3MT\times 1\\ \text{window} ~~~~ \mathbf{W} ~~~ 3MT\times MT\\ \text{means} ~~~~ \boldsymbol{\mu}_{\mathbf{q}} ~~~ 3MT\times 1\\ \text{covariance} ~~~~ \boldsymbol{\Sigma}_\mathbf{q} ~~~ 3MT\times 3MT \end{aligned} 在开始推导参数生成算法之前，我们需要给出线性代数中的一些定理。 \begin{aligned} \frac{dx^{\top}Ax}{dx} = x^{\top}(A+A^{\top})\\ \overset{\text{if}A^{\top}=A}{=}2x^{\top}A^{\top}=2(Ax)^{\top}\\ \frac{dAx}{dx} = A \end{aligned} 二、极大似然参数生成算法（MLPG） 给定高斯分布序列$\mathbf{q}$，参数生成的准则表示为 \begin{aligned} \mathbf{o}_{\max} = \arg\max_{\mathbf{o}}p(\mathbf{o}\mid\mathbf{q},\lambda)\\ = \arg\max_{\mathbf{o}}\mathcal{N}(\mathbf{o}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}}) \end{aligned} 在这里，生成的参数序列既是高斯分布的均值矢量序列，生成的语音变化不自然。为避免这个问题，引入动态窗系数作为约束(Tokuda et al., 2000)。在这个约束条件下，以$\mathbf{o}$为变量的函数最大化就等价于以$\mathbf{c}$为变量的函数最大化。 $$\mathbf{c}_{max}=\arg\max_{\mathbf{c}}\mathcal{N}(\mathbf{Wc}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}})$$ 根据多变量高斯密度函数的定义，得到 \begin{aligned} \mathcal{N}(\mathbf{Wc}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}}) = \frac{1}{\sqrt{(2\pi)^{3MT}|\boldsymbol{\Sigma}_{\mathbf{q}}|}} \exp\left\{-\frac{1}{2}\left(\mathbf{Wc}-\boldsymbol{\mu}_{\mathbf{q}}\right)^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\left(\mathbf{Wc}-\boldsymbol{\mu}_{\mathbf{q}}\right)\right\}\\ = \frac{1}{\sqrt{(2\pi)^{3MT}|\boldsymbol{\Sigma}_{\mathbf{q}}|}} \exp\left\{-\frac{1}{2}\left(\mathbf{c}^{\top}\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}-\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\right)\left(\mathbf{Wc}-\boldsymbol{\mu}_{\mathbf{q}}\right)\right\}\\ = \frac{1}{\sqrt{(2\pi)^{3MT}|\boldsymbol{\Sigma}_{\mathbf{q}}|}} \exp\left\{-\frac{1}{2}\left(\mathbf{c}^{\top}\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}+\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}}-\mathbf{c}^{\top}\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}}-\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)\right\}\\ \left(\because\mathbf{c}^{\top}\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}}=\left(\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)^{\top}=\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)\\ = \frac{1}{\sqrt{(2\pi)^{3MT}|\boldsymbol{\Sigma}_{\mathbf{q}}|}} \exp\left\{-\frac{1}{2}\left(\mathbf{c}^{\top}\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}+\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}}-2\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)\right\} \end{aligned} 似然度定义为密度函数的自然对数。对似然函数求$\mathbf{c}$的偏导数得到 \begin{aligned} \frac{\partial\log\mathcal{N}(\mathbf{Wc}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}})}{\partial\mathbf{c}} = -\frac{1}{2}\left\{\frac{\partial}{\partial\mathbf{c}}\left(\mathbf{c}^{\top}\underline{\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}}\mathbf{c}\right)-2\frac{\partial}{\partial\mathbf{c}}\left(\underline{\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}}\mathbf{c}\right)\right\}\\ \left(\because\left(\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}\right)^{\top}=\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}\right)\\ = -\frac{1}{2}\left\{2\left(\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)^{\top}-2\left(\boldsymbol{\mu}_{\mathbf{q}}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}\right)\right\}\\ = -\frac{1}{2}\left\{2\left(\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{Wc}\right)^{\top}-2\left(\mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}}\right)^{\top}\right\}. \end{aligned} 令$\log\mathcal{N}(\mathbf{Wc}_{\max}\mid\boldsymbol{\mu}_{\mathbf{q}},\boldsymbol{\Sigma}_{\mathbf{q}})/\partial\mathbf{c}_{\max}=\boldsymbol{0}$, 可以得到如下线性方程（组） $$\mathbf{R}_{\mathbf{q}}\mathbf{c}_{\max} = \mathbf{r}_{\mathbf{q}}$$ 其中 \begin{aligned} \mathbf{R}_{\mathbf{q}} = \mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\mathbf{W}\\ \mathbf{r}_{\mathbf{q}} = \mathbf{W}^{\top}\boldsymbol{\Sigma}_{\mathbf{q}}^{-1}\boldsymbol{\mu}_{\mathbf{q}} \end{aligned} 通过求解方程（组）即可得到极大似然准则下的参数序列。 三、实现方法 目前，有三种官方发布的MLPG算法实现， HTS工具包中的HMGenS工具 SPTK工具包中的mlpg工具 hts_engine_API库中的参数生成模块 后两者的应用场景与本文相同，即给定高斯分布序列。在第一个工具中，实现了三种不同的参数生成模式，包括联合最优高斯分布序列的搜索(Tokuda et al., 2000)。但是，通过设置$M=1$，根据HMM的维度间独立假设，三种工具都实现了按维求取。使用Cholesky 分解，$\mathbf{R}_{\mathbf{q}}$可以表示为 $$\mathbf{R}_{\mathbf{q}} = \mathbf{U}_{\mathbf{q}}^{\top}\mathbf{U}_{\mathbf{q}},$$ 其中$\mathbf{U}_{\mathbf{q}}$是一个上三角矩阵。因此参数求解方程（组）可以分解为两个方程（组） \begin{aligned} \mathbf{U}_{\mathbf{q}}^{\top}\mathbf{g}_{\mathbf{q}} = \mathbf{r}_{\mathbf{q}},\\ \mathbf{U}_{\mathbf{q}}\mathbf{c}_{\max} = \mathbf{g}_{\mathbf{q}}. \end{aligned} 上述方程组可以通过前向-后向迭代法（在线性代数课程中，也称为高斯消去法）求解。在以上任一工具的源代码中你都可以找到严格按照上述算法实现的参数求解模块。 目前为止，以上工具包公开发布的稳定版本为 HTS-2.2 SPTK-3.6 hts_engine_API-1.06 参考文献 Keiichi Tokuda, Takao Kobayashi, and Satoshi Imai. speech parameter generation from HMM using dynamic features. In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, pages 660–663, September 1995. Keiichi Tokuda, Takayoshi Yoshimura, Takashi Masuko, Takao Kobayashi, and Tadashi Kitamura. speech parameter generation algorithms for HMM-based speech synthesis. In Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, pages 1315–1318, June 2000. Heiga Zen, Keiichi Tokuda, and Alan W Black. Statistical parametric speech synthesis. Speech Communication, 51(11):1039–1064, November 2009.
个人分类: 课题积累|8355 次阅读|0 个评论
HMM模型的初步学习
friendpine 2013-1-1 11:09
这几天学习了HMM模型，得到了下面的一些粗浅认识，以及对于一个简单的HMM模型的实现（Perl语言实现的）。code写得很冗余，要是把计算过程写成函数会更好。 本文只是自己的一个简单的总结，需要了解的更详细可以参考Lawrence.R. Rabiner的经典论文"A tutorial on Hidden Markov Models and selected applications in Speech recognition"。 1 什么是HMM模型？ 隐马尔科夫模型(Hidden Markov Model),它是一种统计学模型，用来描述随机事件的发生。该模型与马尔科夫模型的差别在于它的状态是隐含的，我们看到的是事件的发生。隐状态之间有转换概率，从隐状态到事件之间有映射概率。 2 HMM模型能够干什么？ 三个问题： 1）给出HMM模型和观测序列，计算在该模型下得到该观测序列的概率 2）给出HMM模型和观测序列，计算在该模型下得到该观测序列最有可能的隐状态序列 3）给定观测序列，估计HMM模型 3 对于每个问题的解答以及例子 对于问题1，一般采用的是forward algorithm(也可以用backward algorithm)。该算法的实现如下（Perl语言）： #指定HMM模型的参数 my $stateNum=3; my $obsNum=10; my @initialP=qw(0.1 0.5 0.4); my @state=qw(A B C); my @obsKind=qw(1 2 3 4); my $obsKind2num; for(my $i=0;$iscalar @obsKind;$i++){ $obsKind2num-{$obsKind }=$i; } my @Atransition=(0.8,0.1,0.1); my @Btransition=(0.2,0.6,0.2); my @Ctransition=(0.3,0.1,0.6); my $transitionRef-{'A'}=\@Atransition; $transitionRef-{'B'}=\@Btransition; $transitionRef-{'C'}=\@Ctransition; my @Atransmit=(0.1,0.2,0.3,0.4); my @Btransmit=(0.2,0.1,0.2,0.5); my @Ctransmit=(0.3,0.4,0.1,0.2); my $transmitRef-{'A'}=\@Atransmit; $transmitRef-{'B'}=\@Btransmit; $transmitRef-{'C'}=\@Ctransmit; my @obs=qw(1 4 2 3 1 4 1 1 4 3); my $proRef; #前向向量 for(my $i=0;$iscalar @obs;$i++){ my $obsKindNum=$obsKind2num-{$obs }; if($i ==0){ for(my $i=0;$iscalar @state;$i++){ $proRef-{0}-{$state }=$initialP *$transmitRef-{$state }- ; } next; } for(my $k=0;$kscalar @state;$k++){ my $pro=0; for(my $h=0;$hscalar @state;$h++){ $pro+=$proRef-{$i-1}-{$state }*$transitionRef-{$state }- ; } $proRef-{$i}-{$state }=$pro*$transmitRef-{$state }- ; } } #把T时刻各个状态的概率相加得到观测序列的概率 my $proFinal=0; foreach my$state(@state){ $proFinal+=$proRef-{(scalar @obs-1)}-{$state}; } print "The final probability is: ",$proFinal,"\n"; 对于问题2,同样是forward algorithm和backword algorithm都可以用。forward algorithm的实现如下： #指定HMM模型的参数 my $stateNum=3; my $obsNum=10; my @initialP=qw(0.1 0.5 0.4); my @state=qw(A B C); my @obsKind=qw(1 2 3 4); my $obsKind2num; for(my $i=0;$iscalar @obsKind;$i++){ $obsKind2num-{$obsKind }=$i; } my @Atransition=(0.8,0.1,0.1); my @Btransition=(0.2,0.6,0.2); my @Ctransition=(0.3,0.1,0.6); my $transitionRef-{'A'}=\@Atransition; $transitionRef-{'B'}=\@Btransition; $transitionRef-{'C'}=\@Ctransition; my @Atransmit=(0.1,0.2,0.3,0.4); my @Btransmit=(0.2,0.1,0.2,0.5); my @Ctransmit=(0.3,0.4,0.1,0.2); my $transmitRef-{'A'}=\@Atransmit; $transmitRef-{'B'}=\@Btransmit; $transmitRef-{'C'}=\@Ctransmit; my @obs=qw(1 1 4 2 2 3 3); my $current2preMax; #用来回溯状态序列 my $proRef; #$proRef-{$index}-{$state} 前向向量 for(my $index=0;$index scalar @obs;$index++){ my $obsNum=$obsKind2num-{$obs }; my $currentObs=$obs ; if($index==0){ for(my $i=0;$iscalar @state;$i++){ $proRef-{$index}-{$state }=$initialP *$transmitRef-{$state }- ; } next; } for(my $i=0;$iscalar @state;$i++){ my $subRef; for(my $j=0;$jscalar @state;$j++){ $subRef-{$state }=$proRef-{$index-1}-{$state }*$transitionRef-{$state }- ; } my @previousState=sort {$subRef-{$b}=$subRef-{$a}}keys %{$subRef}; $proRef-{$index}-{$state }=$subRef-{$previousState }*$transmitRef-{$state }- ; $current2preMax-{$index}-{$state }=$previousState ; } } my @index=sort {$b=$a}keys %{$current2preMax}; my @finalStateSort=sort {$proRef-{$index }-{$b} = $proRef-{$index }-{$a}}keys %{$proRef-{$index }}; print "The final step: ",$index ,"\tthe state with the largest probiligy: ",$finalStateSort ,"\n"; #回溯过程 my @stateSeq; my $finalState=$finalStateSort ; push(@stateSeq,$finalState); for(my $i=$index ;$i0;$i--){ my $cur=$stateSeq ; my $next=$current2preMax-{$i}-{$cur}; push(@stateSeq,$next); } @stateSeq=reverse @stateSeq; print join("-",@stateSeq),"\n"; print join("-",@obs),"\n"; 对于问题3，采用forward-backward algorithm，也称为 Baum-Welch算法。它是一种贪婪算法，往往只能够达到局部最优。它的实现如下： ############################################################################## #采用Forward-backward algorithm （Baum-Welch算法）估计HMM模型的参数； ############################################################################## #模型参数的初始化 my $stateNum=3; my $obsNum=10; my @initialP=qw(0.1 0.5 0.4); my @state=qw(A B C); my @obsKind=qw(1 2 3 4); my $obsKind2num; for(my $i=0;$iscalar @obsKind;$i++){ $obsKind2num-{$obsKind }=$i; } my @Atransition=(0.8,0.1,0.1); my @Btransition=(0.2,0.6,0.2); my @Ctransition=(0.3,0.1,0.6); my $transitionRef-{'A'}=\@Atransition; $transitionRef-{'B'}=\@Btransition; $transitionRef-{'C'}=\@Ctransition; my @Atransmit=(0.1,0.2,0.3,0.4); my @Btransmit=(0.2,0.1,0.2,0.5); my @Ctransmit=(0.3,0.4,0.1,0.2); my $transmitRef-{'A'}=\@Atransmit; $transmitRef-{'B'}=\@Btransmit; $transmitRef-{'C'}=\@Ctransmit; my @obs=qw(1 1 4 1 2 1 1 2 3 4 2 2 3 3 4 2); #观测序列 #初始模型下得到观测序列的概率; my $proBefore=0; my $proRef; for(my $i=0;$iscalar @obs;$i++){ my $obsKindNum=$obsKind2num-{$obs }; if($i ==0){ for(my $i=0;$iscalar @state;$i++){ $proRef-{0}-{$state }=$initialP *$transmitRef-{$state }- ; } next; } for(my $k=0;$kscalar @state;$k++){ my $pro=0; for(my $h=0;$hscalar @state;$h++){ $pro+=$proRef-{$i-1}-{$state }*$transitionRef-{$state }- ; } $proRef-{$i}-{$state }=$pro*$transmitRef-{$state }- ; } } foreach my$state(keys %{$proRef-{(scalar @obs)-1}}){ $proBefore+=$proRef-{(scalar @obs)-1}-{$state}; } $proBefore=log($proBefore); open(OUT,"test")or die"$!"; my $proAfter=$proBefore; my $threshold=1e-10; #判断迭代过程终止的cutoff for(my $time=1;$time 1000;$time++){ my $epsonRef; #$epsonRef-{$t}-{$i}-{$j} my $alphaRef; #$alphaRef-{$t}-{$i} my $betaRef; #$betaRef-{$t}-{$i} #calculate $alphaRef; for(my $t=0;$t scalar @obs;$t++){ my $obsKindNum=$obsKind2num-{$obs }; if($t==0){ for(my $i=0;$iscalar @state;$i++){ $alphaRef-{0}-{$state }=$initialP *$transmitRef-{$state }- ; } next; } for(my $i=0;$iscalar @state;$i++){ my $pro=0; for(my $j=0;$jscalar @state;$j++){ $pro+=$alphaRef-{$t-1}-{$state }*$transitionRef-{$state }- ; } $alphaRef-{$t}-{$state }=$pro*$transmitRef-{$state }- ; } } #calculate $betaRef; for(my $t=(scalar @obs-1);$t =0;$t--){ if($t==((scalar @obs)-1)){ for(my $i=0;$iscalar @state;$i++){ $betaRef-{$t}-{$state }=1; } next; } my $obsKindNum=$obsKind2num-{$obs }; for(my $i=0;$iscalar @state;$i++){ my $pro=0; for(my $j=0;$jscalar @state;$j++){ $pro+=$betaRef-{$t+1}-{$state }*$transitionRef-{$state }- *$transmitRef-{$state }- ; } $betaRef-{$t}-{$state }=$pro; } } #calculate $epsonRef; for(my $t=0;$t scalar @obs-1;$t++){ my $obsKindNum=$obsKind2num-{$obs }; my $denominator=0; for(my $i=0;$iscalar @state;$i++){ for(my $j=0;$jscalar @state;$j++){ $denominator+=$alphaRef-{$t}-{$state }*$transitionRef-{$state }- *$transmitRef-{$state }- *$betaRef-{$t+1}-{$state }; } } for(my $i=0;$iscalar @state;$i++){ for(my $j=0;$jscalar @state;$j++){ my $numerator=0; $numerator+=$alphaRef-{$t}-{$state }*$transitionRef-{$state }- *$transmitRef-{$state }- *$betaRef-{$t+1}-{$state }; $epsonRef-{$t}-{$state }-{$state }=$numerator/$denominator; } } } #calculate $rRef; my $rRef; for(my $t=0;$tscalar @obs;$t++){ for(my $i=0;$iscalar @state;$i++){ my $r=0; foreach my$state(keys %{$epsonRef-{$t}-{$state }}){ $r+=$epsonRef-{$t}-{$state }-{$state}; } $rRef-{$t}-{$state }=$r; } } #Re-estimation of the parameters; for(my $i=0;$iscalar @state;$i++){ $initialP =$rRef-{0}-{$state }; } for(my $i=0;$iscalar @state;$i++){ my $denominator=0; for(my $t=0;$t(scalar @obs)-1;$t++){ $denominator+=$rRef-{$t}-{$state }; } for(my $j=0;$jscalar @state;$j++){ my $numerator=0; for(my $t=0;$t (scalar @obs)-1;$t++){ $numerator+=$epsonRef-{$t}-{$state } -{$state }; } $transitionRef-{$state }- =$numerator/$denominator; } } for(my $i=0;$iscalar @state;$i++){ my $denominator=0; for(my $t=0;$t(scalar @obs);$t++){ $denominator+=$rRef-{$t}-{$state }; } for(my $j=0;$jscalar @obsKind;$j++){ my $numerator=0; for(my $t=0;$tscalar @obs;$t++){ if($obs ==$obsKind ){ $numerator+=$rRef-{$t}-{$state }; } } $transmitRef-{$state }- =$numerator/$denominator; } } #计算更新参数后的HMM模型得到观测序列的概率； my $proRef; $proAfter=0; for(my $i=0;$iscalar @obs;$i++){ my $obsKindNum=$obsKind2num-{$obs }; if($i ==0){ for(my $i=0;$iscalar @state;$i++){ $proRef-{0}-{$state }=$initialP *$transmitRef-{$state }- ; } next; } for(my $k=0;$kscalar @state;$k++){ my $pro=0; for(my $h=0;$hscalar @state;$h++){ $pro+=$proRef-{$i-1}-{$state }*$transitionRef-{$state }- ; } $proRef-{$i}-{$state }=$pro*$transmitRef-{$state }- ; } } foreach my$state(keys %{$proRef-{(scalar @obs)-1}}){ $proAfter+=$proRef-{(scalar @obs)-1}-{$state}; } $proAfter=log($proAfter); my $delta=$proAfter-$proBefore; $proBefore=$proAfter; #输出更新后的HMM参数 print OUT $time,"\t",$delta,"\n"; print OUT "initial: ",join("\t",@initialP),"\n"; print OUT "transition matrix: \n"; print OUT "\t",join("\t",@state),"\n"; for(my $i=0;$iscalar @state;$i++){ print OUT $state ; for(my $j=0;$jscalar @state;$j++){ print OUT "\t",$transitionRef-{$state }- ; } print OUT "\n"; } print OUT "transmit matrix: \n"; print OUT "\t",join("\t",@obsKind),"\n"; for(my $i=0;$iscalar @state;$i++){ print OUT $state ; for(my $j=0;$jscalar @obsKind;$j++){ print OUT "\t",$transmitRef-{$state }- ; } print OUT "\n"; } print OUT "\n"; if($delta $threshold){ print OUT $time,"\t",$delta,"\tOver\n"; last; } } close OUT or die"$!"; 更多语言版本的实现见网友的总结： http://hi.baidu.com/shinchen2/item/86fdadcd67d63c7189ad9ee3 。此外，我在学习的过程中参考了下面这位网友的学习经验( http://www.52nlp.cn/hmm-learn-best-practices-seven-forward-backward-algorithm-5). 他写得很通俗易懂，推荐大家看看。 另外，在R语言中专门有一个包HMM用来实现该模型。
个人分类: 生物信息学与计算生物学|9993 次阅读|0 个评论
[转载]隐藏Markov模型
grating 2012-7-12 08:14
不知道在这里如何编辑公式，只好把我原来的文章做一个引用。 http://jacobyuan.blog.sohu.com/160394296.html
个人分类: 翻译|2428 次阅读|0 个评论
生物信息之ME, HMM, MEMM, CRF
热度 1 hjanime 2012-5-22 14:50
做高端的生物信息理论离不开各种modeling 于是乎漫长的digest之路开始... 最大熵模型 Maximum Entropy 现从一个简单例子看起： 比如华盛顿和维吉利亚都可以作人名和地名，而从语料中只知道p(人名)＝0.6，那么p(华盛顿＝人名)的概率为多少比较好呢？一个直观的想法就是p(华盛顿=人名)=0.3。为什么呢？这就是在满足已有证据的情况下不做任何其他假设，也就是熵最大，这就是最大熵模型的原理。 现在来看模型的定义： 首先，明确模型的目标：给定一个上下文x，估计p(y|x) 接着，从训练样本中我们可以得到一串标注过的样本(x_i, y_i)，其中x_i为上下文，y_i \in Y为类别 然后构造特征函数 f(x,y) = 1 如果x,y满足一些条件，比如x=记者*,y＝人名 0 otherwise 注意x是一个上下文，是个向量，而不是单个词 (最大熵模型里的特征概念不同于模式识别里的特征，这里的特征即特征函数，通常是二值函数，也有直接定义成实数的，比如 jeon-sigir06里直接把f定义为KDE距离，不是很明白那样定义的好处。) 于是模型的约束就是对于所有的特征函数模型中的期望等于样本的期望，即 E_p(f) = E_{\tilde p}(f) 其中 E_p(f) = \sum_{x, y}p(x, y)f(x, y) = \sum_{x, y}p(x)p(y|x)f(x,y) \approx \sum_{x, y} \tilde p(x)p(y|x)f(x,y) \tilde p(f) = \sum_{x, y} \tilde p(x, y)f(x, y), 并且对于任意的x：\sum_y p(y|x) = 1 而模型的熵为 H(p)=-\sum_{x,y} \tilde p(x) p(y|x) \log p(y|x) 在满足约束的情况下，使熵最大，于是问题即求 p* =\argmax_{p \in P} -\sum{x, y} p(y|x)\tilde p(x) \log p(y|x) where P={p(y|x) | \all f_i : \sum_{x,y}p(y|x)\tilde p(x)f_i(x,y) = \sum_{x,y}\tilde p(x,y)f_i(x,y), \all x : \sum_y p(y|x) = 1} 可以证明，模型的最优解的形式为 p(y|x) = exp(\sum_i \lambda_i f_i(x,y)) / Zx where Zx = \sum_y exp(\sum_i \lambda_i f_i(x,y)) 隐马尔可夫模型 Hidden Markov Model 马尔可夫模型实际上是个有限状态机，两两状态间有转移概率；隐马尔可夫模型中状态不可见，我们只能看到输出序列，也就是每次状态转移会抛出个观测值；当我们观察到观测序列后，要找到最佳的状态序列。 设O为观测值，x为隐变量，那么模型要找到让P(O)最大的最佳隐藏状态，而 P(O) = \sum_x P(O|X)P(X) 而其中 P(X)=p(x_1)p(x_{2..n}|x_1) =p(x_1)p(x_2|x_1)p(x_{3..n}|x_1,x_2) …… 根据x_i只与x_{i-1}相关的假设有 P(X)=p(x_1)p(x_2|x_1)p(x_3|x_2)…… 而类似的 P(O|X)=p(o_1|x_{1..n})p(o_{2..n}|o_1x_{1..n}) =p(o_1|x_{1..n})p(o_2|o_1x_{1..n})p(o_{3..n}|o_{1,2},x_{1..n}) …… 根据o_i只与x_i有关的假设有 P(O|X)=p(o_1|x_1)p(o_2|x_2)…… 合起来就是 P(O)=\sum_x p(x_1)p(x_2|x_1)p(o_1|x_1)p(x_3|x_2)p(o_2|x_2)…… 定义向前变量\alpha_i(t)为t时刻以状态S_i结束时的总概率 \alpha_j(t)=\sum_{i=1}^N \alpha_ip(x_{t}=j|x_{t-1}=i)p(o_t=i|x_t=i) 定义向后变量\beta_i(t)为给定当前状态S_i和t时刻情况下观测序列中剩余部分的概率和 \beta_i(t)=\sum_{j=1}^N \p(x_{t}=j|x_{t+1}=i)p(o_{t}=i|x_{t}=i) \beta_j(t+1) 于是观测序列的概率为 P(O, X_t=i) = \alpha_i(t)\beta_i(t) 最佳状态可以由动态规划得到 模型参数可以由EM算法得到 最大熵隐马 Maximum Entropy Markov Model HMM的缺点是根据观测序列决定状态序列，是用联合模型解决条件问题；另外，几乎不可能枚举所有所有可能的观测序列。 而MEMM解决了这些问题。 首先，MEMM和MM或HMM都有本质不同，MEMM估计的是P(S|O)，而MM估计的是P(S)，HMM估计的都是P(O)。 P(S|O)=P(s_1|O)P(s_{2..n}|s_1,O) =P(s_1|O)P(s_2|s_1,O)P(s_{3..n}|s_1,s_2,O) …… 然后根据假设有 P(S|O)=P(s_1|O)P(s_{2..n}|s_1,O) =P(s_1|o_1)P(s_2|s_1,o_2)P(s_{3..n}|s_1,s_2,o_3) …… 重新定义特征函数： a=b,r b是指示函数用于指示当前观测 r是状态值 f_a(o_t, S_t) = 1 if b(o_t) is true and s_t = r 于是约束变为 E_a = \sum_{k=1}^m_{s'}\sum_{s \in S}P(s|s', o_k)f_a(o_k, s) / m_s' = \sum_{k=1}^m_{s'} f_a(o_k, s_k) = F_a 这个目标函数和ME的目标函数实质是一样的 于是解的形式为 P(s|s', o)=exp(\sum_a \lambda_a f_a(o, s)) / Z(o, s') 然后依然采用HMM中的前向后向变量，寻找最佳序列 而实际上得到的序列是由计算 P(s|o) = P(s_0)P(s_1|s_0,o_0)P(s_2|s_1, o_1)…… 得到 条件随机场 Conditional Random Fields MEMM其实是用局部信息去优化全局，会有label bias的问题。比如rib和rob，有如下的状态设计： /r 1i - 2 \ 0 b 3 \r 4o - 5 / 如果训练集合里的ri多于ro，那么ro其实就被无视了…… 所以要用全局优化全局，于是引入CRF，其实CRF就是ME扩展版，并不需要由MRF推出 p(y|x)\propto exp(\sum_i\lumbda_k f_k(y_{i-1}, y_i, x)+\sum_k \lumbda_kg_k(x_k, x)) 其实这个定义并保持MRF性质：MRF中clique于clique是独立的 从这点意义上来看，Lafferty也满水的= = 虽然X，Y在图上不需要有相同的结构，甚至X都不需要有图的结构，目前通常G只是一条链G=(V={1,2, ..., m})，E={(i,i+1)}。 如果G是一棵树，它的团就是每条边和其上的点，这时候 p_\theta(y|x) = exp(\sun_{e\in E, k}\lambda_k f_k(e,y|_e, x)+\sum_{v \in V,k}\mu_k g_k(v, y|_v, x)) x是上下文，y是标注序列，y|_s是y中与子图S相连的部分 如果是最简单的first-order情况，g_k就相当于ME中的特征函数，f_k类似于MEMM中的特征函数，就是说用g_k来指示状态，用f_k来指示状态转移 从优化角度来说，这两类函数是一样的，可以简写为f_j(y_{i-1},y_i,x,i)，于是训练类似ME 当CRF是链状时，一个上下文x的标注可以利用矩阵高效计算 对于长度为n的y，定义n+1个|Y|*|Y|矩阵{M_i(x)|i=1..n+1}，其中 Y是类别个数 M_i(y', y|x) = exp(\sum_j \lambda_j f_j(y', y, x, i)) 这个就是第i个矩阵i,j上的元素，表示在x下从y'转移到y的概率 于是有p(y|x, \lambda)=\multi_{i=1}^{n+1}M_i(y_{i-1},y_i|x) / Zx Zx = \multi_{i=1}^{n+1}M_i(x) Zx只是保证概率之和为1 原来已经有人开始做复杂的情况了……04年cvpr有篇用CRF做图像标注的。 看来易想到的idea一般都会有人做了…… John Lafferty, Andrew McCallum这两个人无比牛啊!! 几乎引领着CRF这一领域的发展：虽然ME不是1拉最先提出来的，但是1拉从96年就开始研究crf相关的东西， 这种牛就是每年N篇ICML+N篇NIPS……而且1看起来还是满ss的，套套看，kaka……Andrew McCallum从1999年开始做ME，那年与John Lafferty合作在IJCAI上发了篇用ME做文本分类的文章，2003年在ICML上发了篇Dynamic CRF，然后又把CRF应用在CV上，也是每年N篇ICML+N篇NIPS的人物，1虽然不如John Lafferty ss，但也还是很可爱的！发现牛人的圈子其实很小，一个领域领头人物也就那么几个…… 原来HMM和MEMM,CRF估的东西就不一样，所以以前怎么推都推不出……= = 不过最近发现 conditional random fields ( CRF s) 应用于 model RNA sequence–structure relationship ， predict post-translational modification sites (labels) in protein sequences 越来越广泛了 看来新PGM总有优势所在呀》》 先总结到这儿，欢迎大牛们指点，
个人分类: 生物信息|7428 次阅读|1 个评论
隐马尔科夫模型（HMM）
热度 9 xiaohai2008 2012-1-10 10:36
很久很久以前（大约2007年下半年），本人在读博士期间，主要参考以下两篇文献（一本生物信息分析方面的书，一本HMM的教程类的文章）： R. Durbin, S. Eddy, A. Krogh and G. Mitchison, 1998. Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press, Cambridge. Lawrence R. Rabiner, 1989. A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proceedings of the IEEE, Vol. 77, No. 2, pp. 257-286. 整理了一篇《An Introduction to Hidden Markov Model》，主要对离散型的HMM进行了系统介绍，包括前向算法、后向算法、Viterbi算法以及参数估计算法（其实是著名的EM算法），包括实现和理论两个方面。 缘于文中很多语句是直接从相关文献copy得到，故打算只做为自己和周围的几个朋友学习参考之用，从未打算扩散，但却被人传到了 百度文库 上，居然下载阅读还需支付两个“财富值”。其实本人并不知道这篇拙文被传到了百度文库，是Ranjan Rajendran告诉我的（他在搜HMM时，找到了这篇文章，无法下载，就通过email联系我，我这时才知道）。 是可忍，孰不可忍。 因此，本人决定将其贴于此处（本文见 HMM.pdf ），供有需要的朋友参考。不过，其中许多内容并非本人原创，请大家注意甄别，给您的学习工作带来的麻烦之处，敬请谅解。
个人分类: 机器学习|10815 次阅读|18 个评论
更多...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-29 01:15

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部