HMM(隐马尔可夫)笔记

（模板匹配）DTW主要作用：

把两段不同长度的语音在时间轴上进行了对齐。

DTW不足：

DTW本质上是一种模板匹配技术，只能进行简单的匹配，而且其参数简单，无法对语音信号的多样性建模，只适用于特定人的小词汇量的语音识别。

事实上，人说话不光是时变过程，频域分布也在变化，随机性很大。不同的人所同样的话，语音也存在较大差异。因此语音识别实际上是一个复杂的问题。

双重随机过程

人的发音包含双重随机过程。

想说什么不确定：即说话内容，具体包含哪些符号（音素或字词）？
怎么说不确定：同样内容发音的观察值差异很大。说快还是说慢？其他差异还包括大声还是笑声？在什么环境说？等等。快慢不同会导致声学特征序列长短不同，大小声和环境不同也会导致声学特征分布不同。

需解决问题

如何描述双重随机过程？
如何通过特征序列识别出符号序列？

这就是HMM要解决的问题。

马尔可夫链

用于描述随机过程。在一个随机过程中，如果每个时间的发生概率仅依赖于上一个事件，则称该过程为马尔可夫过程。

假设随机序列在任意时刻，可以处在状态{s1,s2,⋅⋅⋅,sNs_1,s_2,···,s_Ns1,s2,⋅⋅⋅,sN}，且已有随机序列{q1,q2,⋅⋅⋅,qt−1,qtq_1,q_2,···,q_{t-1},q_tq1,q2,⋅⋅⋅,qt−1,qt}，则产生新的事件qt+1q_{t+1}qt+1的概率为：
P(qt+1∣qt,qt−1,⋯,q1)=P(qt+1∣qt)P(q_{t+1}|q_t,q_{t-1},\cdots,q_1)=P(q_{t+1}|q_t) P(qt+1∣qt,qt−1,⋯,q1)=P(qt+1∣qt)
换句话说，马尔可夫过程只能基于当前事件，预测下一个事件，而和之前或未来事件无关。事件和事件都是离散的马尔可夫过程，称为马尔可夫链。

例子如下图：

一个简化的天气模型，其中“晴天”、“雨天”和“多云”各表示一种状态（也是唯一的观察值）。每个状态可以自身转移，也可转移到其他两个状态转移弧上有对应的概率，入状态“雨天”转到状态“多云”的概率为0.4，转移到状态“晴天”的概率为0.2，“雨天”自身转移的概率为0.4。所有从同一个状态转移的概率和等于1。

根据马尔可夫链的定义，只根据当前时刻的状态，预测下一个时刻的状态。例如，已知今天是“晴天”，则未来三天出现“多云”、“晴天”、“雨天”的概率计算如下：
P(多云∣晴天)P(晴天∣多云)P(雨天∣晴天)=0.3∗0.5∗0.1=0.015P(多云|晴天)P(晴天|多云)P(雨天|晴天)=0.3*0.5*0.1=0.015 P(多云∣晴天)P(晴天∣多云)P(雨天∣晴天)=0.3∗0.5∗0.1=0.015
马尔可夫的每个状态只有唯一的观察事件，即状态与观察事件之间不存在随机性。以上天气模型就是一个马尔可夫模型，观察事件序列直接对应状态之间的转移序列。

HMM包含隐含状态，隐含状态和观察事件并不是一一对应关系，因此它所描述的问题比马尔可夫模型更复杂。

本质上，HMM描述了双重随机过程。

马尔可夫链：状态转移的随机性
依存于状态的观察事件的随机性

天气温度变化与HMM

我们以一天的气温变化为例来介绍hmm。对空气的温度，以人的感觉为标准，可以比较主管地划分为冷、凉爽、温暖、热、酷热等几类，这里称为5个状态。温度是可以用温度计定量地测定的，但温度的数值与这五个状态却不是一一对应的。

人的感觉是不能准确地量化的，比如不能简单地将20度强行称为凉爽或是温暖，但是20度肯定不是酷热或冷。对于这种温度到感觉之间的对应关系，可以用模糊数学中的隶属度来描述。但是在概率模型中，则是用统计的办法，用概率密度函数来描述。比如，20度对“凉爽”的输出概率可能为0.7，而对“温暖”的输出概率可能为0.3.

温度的变化是随机的，它的数值随着事件的变化本身就是一个随机过程。温度给人的感觉所对应的五种状态随事件的变化则是另一个随机过程。可以注意到，温度的变化有时是稳定的，可以认为它主流与某个状态中；也可能忽冷忽热，这时说它发生了从一个状态到另一个状态的转移。

由于当前状态与前面一时刻所处的状态有关，所以这一随机过程为马尔可夫过程。

在温度变化的随机过程中，可以用温度计观测到温度的数值，而在状态变化的随机过程中，却无法准确地指出某一时刻到底是什么状态，也就是说，状态的变化是隐含的，这就是所谓“隐含”的来历。

总之，HMM的精髓就在于：观察可测，状态隐含。这里观察就是温度，状态就是人的感觉。

语音识别中的HMM

在语音识别中，所谓的观察序列就是通过计算得到的一帧帧语音参数，如MFCC参数。而状态则是在训练阶段事先规定好的不同语音单元。对于汉语普通话来说，语音单元可以是一个完整的音节，也可以是声母或韵母，还可以是更为精细的音素。

HMM的定义

对HMM有个实例的理解了，现在就开始给出HMM的概念例子：

图中包含了6个状态，T个观察值序列O={o1,o2,o3,o4,...,oTo_1,o_2,o_3,o_4,...,o_To1,o2,o3,o4,...,oT}。

状态与状态之间转移的可能性用转移概率表示，如状态1到状态2的概率为a12a_{12}a12,状态1自身转移概率为a11a_{11}a11，每个状态的所有转移概率和要等于1。

符号介绍

为描述双重随机过程，HMM包括如下要素：

数字符号	符号解释
N	模型中的状态数目
M	每个状态可能输出的观察符号的数目
S	状态集合S={s1,s2,s3,…,sN}，与状态个数N相关
Q	状态转移序列Q={q1,q2,q3,…,qT}，与调用次数T相关
aija_{ij}aij	状态转移概率
A={aij}A=\{a_{ij}\}A={aij}	状态转移概率分布
bj(k)b_j(k)bj(k)	基于状态S_j得到第K个观察值的概率
B={bj(k)}B=\{b_j(k)\}B={bj(k)}	观察符号的概率分布
π={πi}π=\{π_i\}π={πi}	初始状态的概率分布
O	观察值序列O={o1,o2,o3,o4,…,oT}
λ=(π,A,B)	把ABπ简化成一个模型

HMM过程

当给定模型λ=(π,A,B)后，就可将该模型看成一个符号生成器(或称信号源)，生成观察值序列O=o1,o2,…,oT。生成过程（也称HMM过程）如下：

初试状态概率分布为π，随机选择一个初试状态q1q_1q1=sis_isi。
令t=1。
基于状态SiS_iSi的符号概率分布为bi(k)b_i(k)bi(k)，随机产生一个输出符号ot=Vko_t=V_kot=Vk
基于状态SiS_iSi的状态转移概率分布为aija_{ij}aij，随机转移至一个新的状态qt+1=Sjq_t+1 = S_jqt+1=Sj
令t=t+1，若t<=T,则返回步骤（3），否则结束过程。

状态转移序列Q(q1,q2,…,qT)和观察值序列O(o1,o2,…,oT)一一对应。

状态转移序列记录的是HMM某个状态的索引，它和观察值序列不是一一对应的。

也就是HMM的某个状态可能会多次出现，如s1,s1,s2,s2,s2,s3,s3,…。即代表第1帧和第2帧观察值的q1和q2均对应HMM的状态s1。

HMM的三个基本问题

模型评估问题：如何求概率P(O|λ)
最佳路径问题：如何求隐含状态序列Q=q1,q2,…qT
模型训练问题：如何求模型参数π,A,B

模型评估问题

当给定模型λ=(π,A,B)以及观察序列O，一种办法是采用穷举法，计算模型λ对观察值序列O的P(O|λ)概率，步骤如下：

穷举法

对长度为T的观察值序列O，找出所有可能产生该观察值序列O的状态转移序列Q=q1j,q2j,qTj(j=1,2,⋅⋅⋅,J)Q=q^j_1,q^j_2,q^j_T(j=1,2,···,J)Q=q1j,q2j,qTj(j=1,2,⋅⋅⋅,J)
分别计算QjQ^jQj与观察值序列O的联合概率P(O,Qj∣λ)P(O,Q^j|λ)P(O,Qj∣λ)
取各联合概率P(O,Qj∣λ)P(O,Q^j|λ)P(O,Qj∣λ)的和，即：

P(O∣λ)=∑j=1JP(O,Qj∣λ)P(O|\lambda)=\sum^J_{j=1}P(O,Q^j|\lambda) P(O∣λ)=j=1∑JP(O,Qj∣λ)

将P(O,Q^j|λ)进一步表示为
P(O,Qj∣λ)=P(Qj∣λ)P(O∣Qj,λ)P(O,Q^j|\lambda)=P(Q^j|\lambda)P(O|Q^j,\lambda) P(O,Qj∣λ)=P(Qj∣λ)P(O∣Qj,λ)
分别计算右边两项：
P(Qj∣λ)=p(q1j)p(q2j∣q1j)p(q3j∣q2j)⋯p(qTj∣qT−1j)=a0,1ja1,2ja2,3j⋯aT−1,TjP(Q^j|\lambda)=p(q^j_1)p(q^j_2|q^j_1)p(q^j_3|q^j_2)\cdots p(q^j_T|q^j_{T-1})=a^j_{0,1}a^j_{1,2}a^j_{2,3}\cdots a^j_{T-1,T} P(Qj∣λ)=p(q1j)p(q2j∣q1j)p(q3j∣q2j)⋯p(qTj∣qT−1j)=a0,1ja1,2ja2,3j⋯aT−1,Tj

P(O∣Qj,λ)=p(o1∣q1j)p(o2∣q2j)⋯p(oT∣qTj)=b1j(o1)b2j(o2)b3j(o3)⋯bTj(oT)P(O|Q^j,\lambda)=p(o_1|q^j_1)p(o_2|q^j_2)\cdots p(o_T|q^j_T)=b^j_1(o_1)b^j_2(o_2)b^j_3(o_3)\cdots b^j_T(o_T) P(O∣Qj,λ)=p(o1∣q1j)p(o2∣q2j)⋯p(oT∣qTj)=b1j(o1)b2j(o2)b3j(o3)⋯bTj(oT)

最后得到
P(O,Qj∣λ)=a0,1jb1j(o1)a1,2jb2j(o2)⋯aT−1,TjbTj(oT)P(O,Q^j|\lambda)=a^j_{0,1}b^j_1(o_1)a^j_{1,2}b^j_2(o_2)\cdots a^j_{T-1,T}b^j_T(o_T) P(O,Qj∣λ)=a0,1jb1j(o1)a1,2jb2j(o2)⋯aT−1,TjbTj(oT)

P(O∣λ)=∑j=1JP(O,Qj∣λ)=∑j=1J∏t=1Tat−1,tjbtj(ot)P(O|\lambda)=\sum^J_{j=1}P(O,Q^j|\lambda) =\sum^J_{j=1} \prod^T_{t=1} a^j_{t-1,t}b^j_t(o_t) P(O∣λ)=j=1∑JP(O,Qj∣λ)=j=1∑Jt=1∏Tat−1,tjbtj(ot)

穷举法算起来就比较的慢，时间复杂度高，因此我们还一种更高效的算法。

前向后向算法

前向后向算法用来解决高效计算P(O|λ)的问题。主要分前向和后向两部分。

前向算法：

前向算法按输出观察值序列的时间，从前向后递推计算输出概率。此算法用αt(j)α_t(j)αt(j)表示已经输出的观察值o1,o2,…,ot，并且到达状态sjs_jsj的概率为
αt(j)=P(o1,o2,⋯,ot,qt=sj∣λ)\alpha_t(j)=P(o_1,o_2,\cdots,o_t,q_t=s_j|\lambda) αt(j)=P(o1,o2,⋯,ot,qt=sj∣λ)
步骤：

初始化
α1(i)=πibi(o1),1≤i≤N\alpha_1(i)=\pi_ib_i(o_1),1\leq i\leq N α1(i)=πibi(o1),1≤i≤N
迭代计算
αt+1(j)=[∑i=1Nαt(i)aij]bj(ot+1),1≤t≤T−1,1≤j≤N\alpha_{t+1}(j)=[\sum^N_{i=1}\alpha_t(i)a_{ij}]b_j(o_{t+1}),1\leq t\leq T-1,1\leq j\leq N αt+1(j)=[i=1∑Nαt(i)aij]bj(ot+1),1≤t≤T−1,1≤j≤N
终止计算

P(O∣λ)=∑i=1NαT(i)P(O|\lambda)=\sum_{i=1}^N\alpha_T(i) P(O∣λ)=i=1∑NαT(i)

后向算法：

后向算法由后向前推算输出概率。用βt(i)β_t(i)βt(i)表示时刻t时的状态sis_isi，输出结束时的状态为sNs_NsN,则输出观察值序列为ot,ot+1,...,oTo_t,o_{t+1},...,o_Tot,ot+1,...,oT的概率为：
βt(i)=P(ot,ot+1,⋯,oT,qt=si,qT=sN,λ)\beta_t(i)=P(o_t,o_{t+1},\cdots,o_T,q_t=s_i,q_T=s_N,\lambda) βt(i)=P(ot,ot+1,⋯,oT,qt=si,qT=sN,λ)

初始化
βT(i)=1,1≤i≤N\beta_T(i)=1,1\leq i\leq N βT(i)=1,1≤i≤N
迭代计算

βt(i)=∑j=1Naijbj(ot+1)βt+1(j),1≤t≤T−1,1≤j≤N\beta_t(i)=\sum^N_{j=1}a_{ij}b_j(o_{t+1})\beta_{t+1}(j),1\leq t\leq T-1,1\leq j\leq N βt(i)=j=1∑Naijbj(ot+1)βt+1(j),1≤t≤T−1,1≤j≤N

HMM的前向、后向概率估计

由于一个状态需要联系上下文，可用αt(i)α_t(i)αt(i)和βt(i)β_t(i)βt(i)组合来计算P(O|λ)，这样计算的好处是能够把不同时刻的中间结果保存下来，避免不必要的重复计算。
P(O∣λ)=∑i=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j),1≤t≤T−1P(O|\lambda)=\sum^N_{i=1}\sum^N_{j=1}\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j),1\leq t\leq T-1 P(O∣λ)=i=1∑Nj=1∑Nαt(i)aijbj(ot+1)βt+1(j),1≤t≤T−1

P(O∣λ)=∑i=1Nαt(i)βt(i)=∑i=1NαT(i),1≤t≤T−1P(O|\lambda)=\sum^N_{i=1}\alpha_t(i)\beta_{t}(i)=\sum^N_{i=1}\alpha_T(i),1\leq t\leq T-1 P(O∣λ)=i=1∑Nαt(i)βt(i)=i=1∑NαT(i),1≤t≤T−1

实际计算中首先计算出对于每个t和每个状态i的前向概率和后向概率，然后套用上面的公式，计算出该序列对模型的输出概率。这两个公式也称为全概率公式。

最佳路径问题

Viterbi算法

Viterbi算法用于解决如何寻找与给定观察值序列对应的最佳状态序列的问题。

Viterbi算法描述：

定义最佳状态序列Q∗=q1∗,q2∗,...,qT∗，φt(j)Q^*=q^*_1,q^*_2,...,q^*_T，φ_t(j)Q∗=q1∗,q2∗,...,qT∗，φt(j)记录局部最佳状态序列.
定义δt(i)δ_t(i)δt(i)为截止时刻t，依照状态转移序列q1,q2,...,qtq_1,q_2,...,q_tq1,q2,...,qt，产生出观察值o1,o2,...,oto_1,o_2,...,o_to1,o2,...,ot的最大概率，且最终状态为sis_isi。

δt(i)=max⁡q1,q2,⋯,qt−1P(q1,q2,⋯,qt,qt=si∣λ)\delta_t(i)=\max_{q_1,q_2,\cdots,q_{t-1}}P(q_1,q_2,\cdots,q_t,q_t=s_i|\lambda) δt(i)=q1,q2,⋯,qt−1maxP(q1,q2,⋯,qt,qt=si∣λ)

Viterbi算法步骤：

初始化
δ0(1)=1,δ0(j)=0(j≠1)\delta_0(1)=1,\delta_0(j)=0(j\neq1) δ0(1)=1,δ0(j)=0(j=1)

φ1(j)=q1\varphi_1(j)=q_1 φ1(j)=q1
递推
δt(j)=bj(ot)max⁡1≤i≤Nδt−1(i)aij,1≤t≤T;1≤i≤N\delta_t(j)=b_j(o_t)\max_{1\leq i \leq N}\delta_{t-1}(i)a{ij},1\leq t\leq T;1\leq i\leq N δt(j)=bj(ot)1≤i≤Nmaxδt−1(i)aij,1≤t≤T;1≤i≤N

φt(j)=arg⁡max⁡1≤i≤Nδt−1(i)aij\varphi_t(j)=\arg \max_{1\leq i\leq N}\delta_{t-1}(i)a_{ij} φt(j)=arg1≤i≤Nmaxδt−1(i)aij
终止

Pmax(S,O∣λ)=max⁡1≤i≤NδT(i)P_{max}(S,O|\lambda)=\max_{1\leq i\leq N}\delta_T(i) Pmax(S,O∣λ)=1≤i≤NmaxδT(i)

φT(N)=arg⁡max⁡1≤i≤NδT−1(i)aij\varphi_T(N)=\arg\max_{1\leq i\leq N}\delta_{T-1}(i)a_{ij} φT(N)=arg1≤i≤NmaxδT−1(i)aij

算法终止时φ记录的数据便是最佳状态序列Q*。

模型训练问题

模型训练问题课定义为：给定一个观察值序列O=o1,o2,…,oT，确定一个λ=(π,A,B)，使得P(O|λ)最大，用公式表示为：
λˉ=arg⁡max⁡λP(O∣λ)\bar{\lambda}=\arg\max_\lambda P(O|\lambda) λˉ=argλmaxP(O∣λ)
但没有一种方法能直接估计最佳的λ。因此要寻求替代的方法，即根据观察值序列选取初试模型λ1=(π,A1,B1)，然后求得一组新参数λ2=(π,A2,B2)，保证由P(O|λ2)>P(O|λ1)。重复这个过程，逐步改进模型参数，直至P(O|λ)收敛。

基于状态序列Q，由概率公式
P(O,Q∣λ)=πq0∏t=1Taqt−1qtbqt(ot)P(O,Q|\lambda)=\pi_{q_0}\prod^T_{t=1} a_{q_{t-1}q_t}b_{q_t}(o_t) P(O,Q∣λ)=πq0t=1∏Taqt−1qtbqt(ot)
取对数得到
log⁡P(O,Q∣λ)=log⁡πq0+∑t=1Tlog⁡aqt−1qt+∑t=1Tlog⁡bqt(ot)\log P(O,Q|\lambda)=\log\pi_{q_0}+\sum^T_{t=1}\log a_{q_{t-1}q_t}+\sum^T_{t=1}\log b_{q_t}(o_t) logP(O,Q∣λ)=logπq0+t=1∑Tlogaqt−1qt+t=1∑Tlogbqt(ot)
根据Bayes公式和Jensen不等式，经过一系列转化，可定义辅助函数：
Q(λ,λˉ)=∑QP(O,Q∣λ)log⁡P(O,Q∣λ)=∑i=1NP(O,q0=i∣λˉ)log⁡πi=∑i=1N∑t=1TP(O,qt−1=1,qt=j∣λˉ)log⁡aij=∑i=t=1TP(O,qt=i∣λˉ)log⁡bi(ot)Q(\lambda,\bar\lambda)=\sum_QP(O,Q|\lambda)\log P(O,Q|\lambda) \\ =\sum^N_{i=1}P(O,q_0=i|\bar\lambda)\log\pi_i\\ =\sum^N_{i=1}\sum^T_{t=1}P(O,q_{t-1}=1,q_t=j|\bar\lambda)\log a_{ij}\\ =\sum^T_{i=t=1}P(O,q_t=i|\bar\lambda)\log b_i(o_t) Q(λ,λˉ)=Q∑P(O,Q∣λ)logP(O,Q∣λ)=i=1∑NP(O,q0=i∣λˉ)logπi=i=1∑Nt=1∑TP(O,qt−1=1,qt=j∣λˉ)logaij=i=t=1∑TP(O,qt=i∣λˉ)logbi(ot)
由上式可看出，模型参数xi分别似等于π_i,a_ij,b_i(ot)时，均符合如下函数形式：
F(x)=∑icilog⁡xiF(x)=\sum_ic_i\log x_i F(x)=i∑cilogxi
并有条件限制
∑ixi=1\sum_ix_i=1 i∑xi=1
当：
xi=ciΣkckx_i=\frac{c_i}{\Sigma_kc_k} xi=Σkckci
时，该函数可获得全局最优值。

因此我们可得到模型参数的重估计公式：
aij=∑t=1T−1αt(i)aijbj(ot+1)βt+1(j)∑t=1T−1αt(i)βt(i)=∑t=1T−1ξt(i,j)∑t=1T−1γt(i)a_{ij}=\frac{\sum^{T-1}_{t=1}\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum^{T-1}_{t=1}\alpha_t(i)\beta_t(i)}=\frac{\sum^{T-1}_{t=1}\xi_t(i,j)}{\sum^{T-1}_{t=1}\gamma_t(i)} aij=∑t=1T−1αt(i)βt(i)∑t=1T−1αt(i)aijbj(ot+1)βt+1(j)=∑t=1T−1γt(i)∑t=1T−1ξt(i,j)

bi(k)=∑t=1s.t.ot=vkTαt(i)βt(i)∑t=1Tαt(i)βt(i)=∑t=1−s.t.ot=vkTγt(j)∑t=1Tγt(j)b_i(k)=\frac{\sum^{T}_{{t=1}_{s.t.o_t=v_k}}\alpha_t(i)\beta_{t}(i)}{\sum^{T}_{t=1}\alpha_t(i)\beta_t(i)}=\frac{\sum^{T}_{{t=1}-{s.t.o_t=v_k}}\gamma_t(j)}{\sum^{T}_{t=1}\gamma_t(j)} bi(k)=∑t=1Tαt(i)βt(i)∑t=1s.t.ot=vkTαt(i)βt(i)=∑t=1Tγt(j)∑t=1−s.t.ot=vkTγt(j)

可直观地认为aija_{ij}aij的重估计公式是所有时刻从状态sis_isi转移到状态sjs_jsj的概率和除以所有时刻处于状态sis_isi的概率和。

而ξt(i,j)\xi_t(i,j)ξt(i,j)为给定训练序列O和模型λ是，HMM在时刻t处于状态Si，在t+1时刻处于状态sjs_jsj的概率，即：
ξt(i,j)=P(qt=si,qt+1=sj∣O,λ)\xi_t(i,j)=P(q_t=s_i,q_{t+1}=s_j|O,\lambda) ξt(i,j)=P(qt=si,qt+1=sj∣O,λ)
根据前向后向算法，可推导出：
ξt(i,j)=αt(i)aijbj(ot+1)βt+1(j)∑i=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j)=αt(i)aijbj(ot+1)βt+1(j)P(O∣λ)\xi_t(i,j)=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{\sum^{N}_{i=1}\sum^{N}_{j=1}\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}=\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{P(O|\lambda)} ξt(i,j)=∑i=1N∑j=1Nαt(i)aijbj(ot+1)βt+1(j)αt(i)aijbj(ot+1)βt+1(j)=P(O∣λ)αt(i)aijbj(ot+1)βt+1(j)
定义γt(i)\gamma_t(i)γt(i)为在时刻t时处于状态sis_isi的概率：
γt(i)=∑j=1Nξt(i,j)=∑j=1Nαt(i)aijbj(ot+1)βt+1(j)P(O∣λ)=αt(i)βt(i)P(O∣λ)\gamma_t(i)=\sum_{j=1}^N\xi_t(i,j)=\sum_{j=1}^N\frac{\alpha_t(i)a_{ij}b_j(o_{t+1})\beta_{t+1}(j)}{P(O|\lambda)}=\frac{\alpha_t(i)\beta_t(i)}{P(O|\lambda)} γt(i)=j=1∑Nξt(i,j)=j=1∑NP(O∣λ)αt(i)aijbj(ot+1)βt+1(j)=P(O∣λ)αt(i)βt(i)

GMM模型

gmm（高斯混合模型）是采用多个高斯函数组合进行对复杂的数据分布的表示。

K阶高斯GMM使用K个单高斯分布的线性组合来描述。令λ={μ，Σ}，则K阶高斯GMM的概率密度函数如下：
p(x∣λ)=∑k=1Kp(x,k,∣λ)=∑k=1Kp(k)p(x∣k,λ)=∑k=1KckN(x∣μk,Σk)p(x|\lambda)=\sum_{k=1}^Kp(x,k,|\lambda)=\sum^K_{k=1}p(k)p(x|k,\lambda)=\sum^K_{k=1}c_kN(x|\mu_k,\Sigma_k) p(x∣λ)=k=1∑Kp(x,k,∣λ)=k=1∑Kp(k)p(x∣k,λ)=k=1∑KckN(x∣μk,Σk)
c_k是第k个高斯的权重，而∑k+1Kck=1\sum_{k+1}^Kc_k=1∑k+1Kck=1表示所有高斯的权重和为1。第k个高斯函数表示为：
N(x∣μk,Σk)=1(2π)D/2∣Σk∣1/2exp{−(x−μk)TΣk−1(x−μk)2}N(x|\mu_k,\Sigma_k)=\frac{1}{(2\pi)^{D/2}|\Sigma_k|^{1/2}}exp\lbrace-\frac{{(x-\mu_k)^T\Sigma^{-1}_k(x-\mu_k)}}{2}\rbrace N(x∣μk,Σk)=(2π)D/2∣Σk∣1/21exp{−2(x−μk)TΣk−1(x−μk)}

GMM训练步骤

GMM的训练主要分为两步：

初始化，即构造初始模型
重估计，即通过EM迭代算法精细化初始模型

初始化

所有的聚类算法都可用于此，如用K-means、LBG等。

重估计

利用EM的思想，对GMM取对数再对各个参数求偏导置零更新参数数值，根据权重系数ckc_kck限制，加入拉格朗日算子：
ln⁡p(x∣c,μ,Σ)+λ(∑k=1Kck−1)=∑n=1Nln⁡{∑k=1KckN(xn∣μ,Σk)}+λ(∑k=1Kck−1)\ln p(x|c,\mu,\Sigma)+\lambda(\sum^K_{k=1}c_k-1)\\ =\sum_{n=1}^N\ln \{\sum^K_{k=1}c_kN(x_n|\mu,\Sigma_k)\}+\lambda(\sum^K_{k=1}c_k-1) lnp(x∣c,μ,Σ)+λ(k=1∑Kck−1)=n=1∑Nln{k=1∑KckN(xn∣μ,Σk)}+λ(k=1∑Kck−1)

分别对μk、Σk、ck\mu_k、\Sigma_k、c_kμk、Σk、ck求最大似然函数（求导并令导数置0），得：

μk=∑n=1Nγ(n,k)xn∑n=1Nγ(n,k)Σk=∑n=1Nγ(n,k)(xn−μk)(xn−μk)T∑n=1Nγ(n,k)ck=∑n=1Nγ(n,k)∑n=1N∑k=1Kγ(n,k)\mu_k=\frac{\sum_{n=1}^N\gamma(n,k)x_n}{\sum_{n=1}^N\gamma(n,k)}\\ \Sigma_k=\frac{\sum_{n=1}^N\gamma(n,k)(x_n-\mu_k)(x_n-\mu_k)^T}{\sum_{n=1}^N\gamma(n,k)}\\ c_k=\frac{\sum_{n=1}^N\gamma(n,k)}{\sum_{n=1}^N\sum_{k=1}^K\gamma(n,k)}\\ μk=∑n=1Nγ(n,k)∑n=1Nγ(n,k)xnΣk=∑n=1Nγ(n,k)∑n=1Nγ(n,k)(xn−μk)(xn−μk)Tck=∑n=1N∑k=1Kγ(n,k)∑n=1Nγ(n,k)
其中：
γ(n,k)=ckN(xn∣μk,Σk)∑k=1KckN(xn∣μk,Σk)\gamma(n,k)=\frac{c_kN(x_n|\mu_k,\Sigma_k)}{\sum^K_{k=1}c_kN(x_n|\mu_k,\Sigma_k)} γ(n,k)=∑k=1KckN(xn∣μk,Σk)ckN(xn∣μk,Σk)

具体定义gmm步骤：

GMM-HMM统计模型

gmm-hmm模型就是在状态的发射概率b换成了由gmm模型导出到观察值的概率。若观察概率为混合高斯分布形式，即:
bj(ot)=∑k=1KcjkN(ot,μjk,Σjk)b_j(o_t)=\sum^K_{k=1}c_{jk}N(o_t,\mu_jk,\Sigma_{jk}) bj(ot)=k=1∑KcjkN(ot,μjk,Σjk)
模型如下图：

一个状态对应一个GMM，一个GMM有四个高斯。下图为状态2的GMM的参数。

我们也知道，模型训练问题就是通过已有的语音对gmm-hmm模型进行训练，将gmm-hmm进行重估计直至模型收敛即可。

需要重估的参数：

起始概率\pi
转移概率a_ij
各状态中不同pdf的权重c_k
各状态中不同pdf的均值和方差\mu \Sigma

四个参数的估计，均跟前向和后向概率关联，是一种软判决，即分配到状态的观察值多少是用概率来调节的。这样的分配机制跟普通 HMM的训练过程是类似的，即均为Baum-Welch算法。

如果通过Viterbi对齐，是一种硬判决，即每个观察值只跟一个状态对应。

Baum-Welch代码：

function hmm=baum_welch(hmm,obs)mix=hmm.mix; %高斯混合模型
N=hmm.N; %HMM的状态数
K=length(obs); %训练数据样本数
SIZE=size(obs(1).fea,2); %特征矢量的个数for loop = 1:40% ----计算前向, 后向概率矩阵for k=1:Kparam(k)=getparam(hmm,obs(k).fea);end%----重估转移概率矩阵Afor i=1:N-1demon=0;for k=1:K  tmp=param(k).ksai(:,i,:);demon=demon+sum(tmp(:)); %对时间t，j求和endfor j=i:i+1  nom=0;for k=1:K  tmp=param(k).ksai(:,i,j);nom=nom+sum(tmp(:));  %对时间t求和endhmm.trans(i,j)=nom/demon;endend%----重估输出观测值概率Bfor j=1:N %状态循环for l=1:hmm.M(j) %混合高斯的数目%计算各混合成分的均值和协方差矩阵nommean=zeros(1,SIZE);nomvar=zeros(1,SIZE);denom=0;for k=1:K  %训练数目的循环T=size(obs(k).fea,1);  %帧数for t=1:T   %帧数（时间）的遍历x=obs(k).fea(t,:);nommean=nommean+param(k).gama(t,j,l)*x;nomvar=nomvar+param(k).gama(t,j,l)*(x-mix(j).mean(l,:)).^2;denom=denom+param(k).gama(t,j,l);endendhmm.mix(j).mean(l,:)=nommean/denom;hmm.mix(j).var(l,:)=nomvar/denom;%计算各混合成分的权重nom=0;denom=0;for k=1:Ktmp=param(k).gama(:,j,l);nom=nom+sum(tmp(:));tmp=param(k).gama(:,j,:);denom=denom+sum(tmp(:));endhmm.mix(j).weight(l)=nom/denom;endend   end