- Visual C++数字图像模式识别典型案例详解
- 冯伟兴 梁洪 王臣业编著
- 3406字
- 2025-03-16 03:50:24
2.2.2 隐马尔可夫模型
隐马尔可夫模型（Hidden Markov Models，HMM）是序列数据处理和统计学习的一种重要概率模型，具有建模简单、数据计算量小、运行速度快、识别率高等特点，近几年来已经被成功应用到许多工程任务中。
1.隐马尔可夫概念
隐马尔可夫模型（HMM）由两个序列组成，一个是不可观察的（隐藏的）状态变化序列，一个是由该隐藏的状态所产生的可观察的符号序列。隐马尔可夫模型由其初始分布、转移概率矩阵和发射概率矩阵完全确定。一个建好的隐马尔可夫模型的运作过程可以这样描述：由初始分布随机地产生一个状态，然后按转移概率矩阵随机地从一个状态转移到另一个状态，在每一个状态按发射概率矩阵随机地生成一个字符，结果是产生一个状态序列和一个相应的字符序列。
一般情况下，不仅模型能够产生的状态序列不是唯一的，每一个状态序列能够产生的字符序列也不是唯一的。模型以一定的概率产生特定的状态序列，该状态序列又以一定的概率产生特定的字符序列。在实际问题中，通常字符序列是可观测的，状态序列是不可观测的，是“隐藏的”。
HMM可以表示为λ=(A,S,π)，描述如下：
（a）隐状态集S。S={S1,S2,S3,…,SN}，N表示状态数。
（b）观察符号集V。V ={V1,V2,V3,…,VM}，M表示符号数。
（c）状态转换概率矩阵A。A={ai}j
式中
式中，Q=(q1,q2,…,qL)表示状态序列，其中qi∈S，1≤i≤L，L表示序列的长度。t为当前时刻。aij为从状态si转换到s j的概率。
（d）符号释放概率矩阵E。E={e jk}，其中
式中，O=(o1,o2,…,oL)表示符号序列，其中o j∈V，1≤ j≤L。e jk表示在状态S j下产生符号Vk的概率。
（e）初始状态分布π。π={πi}，其中
表示模型起始于状态Si（即时间为0时的状态）的概率。
基于HMM的模式识别方法的基本思想：对于一个可用HMM描述的模式识别问题，其可能出现的所有不同的隐马尔可夫模型λ1,λ2,λ3,…,λn构成一个集合Λ，其中任意一个元素λi为特定的隐马尔可夫模型，其中包括此特定隐马尔可夫模型的转移矩阵、初分布及观测概率。观测概率的定义如下：
给定一组样本观测值(o1,o2,…,on)，选择λ∗∈Λ，使得
则λ∗是观测(o1,o2,…,on)的最可能模式。
2.隐马尔可夫模型基本算法
HMM需解决的基本问题：
·概率问题。给定一个HMM（λ=(A,S,π)）和一个观测符号序列O，求由该模型得到序列O的概率P(O|λ)。
·解码问题：给定HMM（λ=(A,S,π)）和一个观测序列O，解码问题就是寻找最可能的隐状态序列，即如何选择产生该观察序列O的最佳状态序列Q*，也就是如何最好地解释符号序列O。
·学习问题。给定一组观测序列样本集，如何确定其对应的HMM（λ=(A,S,π)），使该模型能够最好地描述观测样本集，也就是该模型的参数能够使P(O|λ)最大。
（1）前向后向算法
对于可能性问题，引入前向变量和后向变量进行动态规划求解。
1）前向算法
前向算法的简单图解如图2-8所示。
图2-8 前向算法示意图
前向变量αj(t)表示在t时刻状态为s j，且序列的前t个符号为o1,o2,…,ot的概率，即
根据前向变量的定义，可以得到：
前向变量的迭代公式：
前向算法的流程：
（a）初始化：
ai(0)=πiei(o 1),1≤i≤N
（b）递归：
t=1,2,3,…,L
j=1,2,3,…,N
a j(t)=ej(ot)∑k(ak(t-1)akj)
next
next
（c）终止：
P (o|λ)=∑kαk(L)
（d）算法结束。
2）后向算法
后向算法的简单图解如图2-9所示。
图2-9 后向算法示意图
后向变量βj(t)表示在t时刻状态为s j，且序列从t+1到L个符号为ot+1,ot+2,…,ot+L的概率，即
根据后向变量的定义，可以得到：
后向变量的迭代公式：
后向算法的流程如下：
（a）初始化：
βi (L)=1 1≤i≤N
（b）递归：
t=L-1,…,1
j= 1,2,…,N
βj(t)=∑kαjk ek(ot+ 1)βk(t+1)
next
next
（c）终止：
Po λ( | )=∑k k k(1) k(1)
（d）算法结束。前向算法和后向算法的时间复杂度为o(N 2L)，空间复杂度为o(NL)。
（2）解码问题
根据给定的HMM(λ=(A,E,π))和符号序列O=(o1,o2,…,oL)，确定一条具有最大概率的状态序列，即最优路径Q*：
应用Viterbi算法求解最优路径。定义变量δj(t)，它表示由状态序列q1,q 2,…,qt且qt =s j产生出符号序列o1,o2,…,ot的最大概率，即
δj (t)可以通过迭代计算得到：
从初始状态开始迭代计算δj(t)，一直到序列末端，然后再回溯得到最优路径。
Viterbi算法流程如下：
（a）初始化：
δi(1)=πiei(o 1) 1≤i≤N
（b）递归：
（c）终止：
（d）回溯：
qt *=φt+1(qt+1)，t=L-1,L-2,…,1
（e）算法结束。
Viterbi算法提供了有效的计算方法来分析HMM中的观察序列，以此得到最可能的隐状态序列，并利用递归来减小计算负担。
（3）训练问题当只有符号序列o1,o2,…,on是给定的时，才要通过训练学习估计模型的各项参数，使观测的符号序列能够“最好”地被描述，即λ能够使P(o1,o2,…,on|λ)最大。
如果符号序列是相互独立的，则它们的联合分布概率即为它们单个概率的乘积：
如果路径已知，则要统计各个特定的初始状态，转移概率以及产生概率。记某种初始状态、转移和出现的次数为πi、Aij和Ei(Vk)，则初始状态概率、转移概率和产生概率的极大似然概率估计值为
如果路径未知，常采用Baum-Welch算法来实现对HMM的参数估计。
根据前向变量和后向变量的定义，可以得到：
通过对所有位置以及所有序列的求和，可以得到：
再根据式（2-116）～式（2-118）即可以得到模型的各项参数估计。
Baum-Welch算法流程如下：
定义变量γt(i)=P(qt =si|O,λ)，ξt(i,j)=P(qt =si,qt+1=s j|O,λ)
则
且
。
故
=时刻t=1在状态Si的期望值=γ1(i)。
此时有
且
，1≤i≤N,1≤ j≤N。
可以通过证明得出结论：
且P(o|λ)迭代收敛到局部最佳。
3.隐马尔可夫模型的C语言实现
下面介绍隐马尔可夫模型的C语言实现代码。
动态数组申请的代码如代码2-2所示。
代码2-2动态数组申请代码
float ＊vector(long nl,long nh)
{
float ＊v;
v=(float ＊)calloc((unsigned) (nh-nl+1),sizeof(float));
if (!v) nrerror("allocation failure in vector()");
return v-nl;
}
int ＊ivector(long nl,long nh)
{
int ＊v;
v=(int ＊)calloc((unsigned) (nh-nl+1),sizeof(int));
if (!v) nrerror("allocation failure in ivector()");
return v-nl;
}
unsigned char ＊cvector(long nl,long nh)
{
unsigned char ＊v;
v=(unsigned char ＊)calloc((unsigned) (nh-nl+1),sizeof(unsigned char));
if (!v) nrerror("allocation failure in cvector()");
return v-nl;
}
unsigned long ＊lvector(long nl,long nh)
{
unsigned long ＊v;
v=(unsigned long ＊)calloc((unsigned) (nh-nl+1),sizeof(unsigned long));
if (!v) nrerror("allocation failure in lvector()");
return v-nl;
}
double ＊dvector(long nl,long nh)
{
double ＊v;
v=(double ＊)calloc((unsigned) (nh-nl+1),sizeof(double));
if (!v) nrerror("allocation failure in dvector()");
return v-nl;
}
float ＊＊matrix(long nrl,long nrh,long ncl,long nch)
{
int i;
float ＊＊m;
m=(float ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(float＊));
if (!m) nrerror("allocation failure 1 in matrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(float ＊) calloc((unsigned) (nch-ncl+1),sizeof(float));
if (!m[i]) nrerror("allocation failure 2 in matrix()");
m[i] -= ncl;
}
return m;
}
double ＊＊dmatrix(long nrl,long nrh,long ncl,long nch)
{
int i;
double ＊＊m;
m=(double ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(double＊));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(double ＊) calloc((unsigned) (nch-ncl+1),sizeof(double));
if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
m[i] -= ncl;
}
return m;
}
int ＊＊imatrix(long nrl,long nrh,long ncl,long nch)
{
int i,＊＊m;
m=(int ＊＊)calloc((unsigned) (nrh-nrl+1),sizeof(int＊));
if (!m) nrerror("allocation failure 1 in imatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(int ＊)calloc((unsigned) (nch-ncl+1),sizeof(int));
if (!m[i]) nrerror("allocation failure 2 in imatrix()");
m[i] -= ncl;
}
return m;
}
unsigned char ＊＊cmatrix(long nrl,long nrh,long ncl,long nch)
{
int i;
unsigned char ＊＊m;
m=(unsigned char ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(unsigned char＊));
if (!m) nrerror("allocation failure 1 in dmatrix()");
m -= nrl;
for(i=nrl;i<=nrh;i++) {
m[i]=(unsigned char ＊) calloc((unsigned) (nch-ncl+1),sizeof(unsigned char));
if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
m[i] -= ncl;
}
return m;
}
//求一个二维矩阵的子矩阵
float ＊＊submatrix(float ＊＊a,long oldrl,long oldrh,long oldcl,long oldch,long
newrl,long newcl)
{
int i,j;
float ＊＊m;
m=(float ＊＊) calloc((unsigned) (oldrh-oldrl+1),sizeof(float＊));
if (!m) nrerror("allocation failure in submatrix()");
m -= newrl;
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;
return m;
}
//将一个实数一维数组转化为二维数组
float ＊＊convert_matrix(float ＊a,long nrl,long nrh,long ncl,long nch)
{
int i,j,nrow,ncol;
float ＊＊m;
nrow=nrh-nrl+1;
ncol=nch-ncl+1;
m = (float ＊＊) calloc((unsigned) (nrow),sizeof(float＊));
if (!m) nrerror("allocation failure in convert_matrix()");
m -= nrl;
for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol＊i-ncl;
return m;
}
前向算法的代码如代码2-3所示。
代码2-3前向算法代码
//前向算法估计参数
void Forward(HMM ＊phmm,int T,int ＊O,double ＊＊alpha,double ＊pprob)
{
int i,j;
int t;
double sum;
//初始化
for(i=1;i<=phmm->N;i++)
alpha[1][i]=phmm->pi[i]＊phmm->B[i][O[1]];
//递归
for(t=1;t<T;t++){
for(j=1;j<=phmm->N;j++){
sum=0.0;
for(i=1;j<=phmm->N;i++)
sum+=alpha[t][i]＊(phmm->A[i][j]);
alpha[t+1][j]=sum＊(phmm->B[j][O[t+1]]);
}
}
//终止
＊pprob=0.0;
for(i=1;i<=phmm->N;i++)
＊pprob+=alpha[T][i];
}
void ForwardWithScale(HMM＊phmm,int T,int ＊O,double ＊＊alpha,double ＊scale,double
＊pprob)
// pprob是对数概率
{
int i,j;
int t;
double sum;
//初始化
scale[1]=0.0;
for(i=1;i<=phmm->N;i++){
alpha[1][i]=phmm->pi[i]＊(phmm->B[i][O[1]+1]);
scale[1]+=alpha[1][i];
}
for(i=1;i<=phmm->N;i++)
alpha[1][i]/=scale[1];
//递归
for(t=1;t<T;t++){
scale[t+1]=0.0;
for(j=1;j<=phmm->N;j++){
sum=0.0;
for(i=1;i<=phmm->N;i++)
sum+=alpha[t][i]＊(phmm->A[i][j]);
alpha[t+1][j]=sum＊(phmm->B[j][O[t+1]]);
scale[t+1]+=alpha[t+1][j];
}
for(j=1;j<=phmm->N;j++)
alpha[t+1][j]/=scale[t+1];
}
//终止
＊pprob=0.0;
for(t=1;t<=T;t++)
＊pprob+=log(scale[t]);
}
后向算法的代码如代码2-4所示。
代码2-4后向算法代码
//后向算法估计参数
void Backward(HMM ＊phmm, int T, int ＊O, double ＊＊beta, double ＊pprob)
{
int i, j;
int t;
double sum;
//初始化
for (i = 1; i <= phmm->N; i++)
beta[T][i] = 1.0;
//递归
for (t = T -1; t >= 1; t--)
{
for (i = 1; i <= phmm->N; i++)
{
sum = 0.0;
for (j = 1; j <= phmm->N; j++)
sum += phmm->A[i][j] ＊ (phmm->B[j][O[t+1]])＊beta[t+1][j];
beta[t][i] = sum;
}
}
//终止
＊pprob = 0.0;
for (i = 1; i <= phmm->N; i++)
＊pprob += beta[1][i];
}
void BackwardWithScale(HMM ＊phmm, int T, int ＊O, double ＊＊beta, double ＊scale,double ＊pprob)
{
int i, j;
int t;
double sum;
//初始化
for (i = 1; i <= phmm->N; i++)
beta[T][i] = 1.0/scale[T];
//递归
for (t = T -1; t >= 1; t--)
{
for (i = 1; i <= phmm->N; i++)
{
sum = 0.0;
for (j = 1; j <= phmm->N; j++)
sum += phmm->A[i][j] ＊ (phmm->B[j][O[t+1]])＊beta[t+1][j];
beta[t][i] = sum/scale[t];
}
}
}
Viterbi算法的代码如代码2-5所示。
代码2-5 Viterbi算法代码
#define VITHUGE 100000000000.0
void Viterbi(HMM ＊phmm, int T, int ＊O, double ＊＊delta, int ＊＊psi, int ＊q, double ＊pprob)
{
int i, j;
int t;
int maxvalind;
double maxval, val;
//初始化
for (i = 1; i <= phmm->N; i++)
{
delta[1][i] = phmm->pi[i] ＊ (phmm->B[i][O[1]]);
psi[1][i] = 0;
}
//递归
for (t = 2; t <= T; t++)
{
for (j = 1; j <= phmm->N; j++)
{
maxval = 0.0;
maxvalind = 1;
for (i = 1; i <= phmm->N; i++)
{
val = delta[t-1][i]＊(phmm->A[i][j]);
if (val > maxval) {
maxval = val;
maxvalind = i;
}
}
delta[t][j] = maxval＊(phmm->B[j][O[t]]);
psi[t][j] = maxvalind;
}
}
//终止
＊pprob = 0.0;
q[T] = 1;
for (i = 1; i <= phmm->N; i++)
{
if (delta[T][i] > ＊pprob)
{
＊pprob = delta[T][i];
q[T] = i;
}
}
for (t = T -1; t >= 1; t--)
q[t] = psi[t+1][q[t+1]];
}
void ViterbiLog(HMM ＊phmm, int T, int ＊O, double ＊＊delta, int ＊＊psi, int ＊q, double
＊pprob)
{
int i, j;
int t;
int maxvalind;
double maxval, val;
double ＊＊biot;
//预处理
for (i = 1; i <= phmm->N; i++)
phmm->pi[i] = log(phmm->pi[i]);
for (i = 1; i <= phmm->N; i++)
for (j = 1; j <= phmm->N; j++)
{
phmm->A[i][j] = log(phmm->A[i][j]);
}
biot = dmatrix(1, phmm->N, 1, T);
for (i = 1; i <= phmm->N; i++)
for (t = 1; t <= T; t++)
{
biot[i][t] = log(phmm->B[i][O[t]]);
}
//初始化
for (i = 1; i <= phmm->N; i++)
{
delta[1][i] = phmm->pi[i] + biot[i][1];
psi[1][i] = 0;
}
//递归
for (t = 2; t <= T; t++)
{
for (j = 1; j <= phmm->N; j++)
{
maxval = -VITHUGE;
maxvalind = 1;
for (i = 1; i <= phmm->N; i++)
{
val = delta[t-1][i] + (phmm->A[i][j]);
if (val > maxval)
{
maxval = val;
maxvalind = i;
}
}
delta[t][j] = maxval + biot[j][t];
psi[t][j] = maxvalind;
}
}
//终止
＊pprob = -VITHUGE;
q[T] = 1;
for (i = 1; i <= phmm->N; i++)
{
if (delta[T][i] > ＊pprob)
{
＊pprob = delta[T][i];
q[T] = i;
}
}
//回朔
for (t = T -1; t >= 1; t--)
q[t] = psi[t+1][q[t+1]];
}
Baum-Welch算法的代码如代码2-6所示。
代码2-6 Baum-Welch算法代码
#define DELTA 0.001
void BaumWelch(HMM ＊phmm, int T, int ＊O, double ＊＊alpha, double ＊＊beta, double＊＊gamma, int ＊pniter, double ＊plogprobinit, double ＊plogprobfinal)
{
int i, j, k;
int t, l = 0;
double logprobf, logprobb;
double numeratorA, denominatorA;
double numeratorB, denominatorB;
double ＊＊＊xi, ＊scale;
double delta, deltaprev, logprobprev;
deltaprev = 10e-70;
xi = AllocXi(T, phmm->N);
scale = dvector(1, T);
ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
＊plogprobinit = logprobf;
BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
ComputeGamma(phmm, T, alpha, beta, gamma);
ComputeXi(phmm, T, O, alpha, beta, xi);
logprobprev = logprobf;
do{
//重新估计t=1时，状态为i的频率
for (i = 1; i <= phmm->N; i++)
phmm->pi[i] = .001 + .999＊gamma[1][i];
//重新估计转移矩阵和观察矩阵
for (i = 1; i <= phmm->N; i++)
{
denominatorA = 0.0;
for (t = 1; t <= T -1; t++)
denominatorA += gamma[t][i];
for (j = 1; j <= phmm->N; j++)
{
numeratorA = 0.0;
for (t = 1; t <= T -1; t++)
numeratorA += xi[t][i][j];
phmm->A[i][j] = .001 + .999＊numeratorA/denominatorA;
}
denominatorB = denominatorA + gamma[T][i];
for (k = 1; k <= phmm->M; k++)
{
numeratorB = 0.0;
for (t = 1; t <= T; t++)
{
if (O[t] == k)
numeratorB += gamma[t][i];
}
phmm->B[i][k] = .001 +.999＊numeratorB/denominatorB;
}
}
ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
ComputeGamma(phmm, T, alpha, beta, gamma);
ComputeXi(phmm, T, O, alpha, beta, xi);
//计算两次直接的概率差
delta = logprobf - logprobprev;
logprobprev = logprobf;
l++;
}
while (delta > DELTA);
＊pniter = l;
＊plogprobfinal = logprobf; /＊ log P(O|estimated model) ＊/
FreeXi(xi, T, phmm->N);
free_dvector(scale, 1, T);
}
void ComputeGamma(HMM ＊phmm, int T, double ＊＊alpha, double ＊＊beta, double ＊＊gamma)
{
int i, j;
int t;
double denominator;
for (t = 1; t <= T; t++)
{
denominator = 0.0;
for (j = 1; j <= phmm->N; j++)
{
gamma[t][j] = alpha[t][j]＊beta[t][j];
denominator += gamma[t][j];
}
for (i = 1; i <= phmm->N; i++)
gamma[t][i] = gamma[t][i]/denominator;
}
}
void ComputeXi(HMM＊ phmm, int T, int ＊O, double ＊＊alpha, double ＊＊beta, double ＊＊＊xi)
{
int i, j;
int t;
double sum;
for (t = 1; t <= T -1; t++) {
sum = 0.0;
for (i = 1; i <= phmm->N; i++)
for (j = 1; j <= phmm->N; j++)
{
xi[t][i][j] = alpha[t][i]＊beta[t+1][j]
＊(phmm->A[i][j])
＊(phmm->B[j][O[t+1]]);
sum += xi[t][i][j];
}
for (i = 1; i <= phmm->N; i++)
for (j = 1; j <= phmm->N; j++)
xi[t][i][j] /= sum;
}
}
double ＊＊＊ AllocXi(int T, int N)
{
int t;
double ＊＊＊xi;
xi = (double ＊＊＊) malloc(T＊sizeof(double ＊＊));
xi --;
for (t = 1; t <= T; t++)
xi[t] = dmatrix(1, N, 1, N);
return xi;
}
void FreeXi(double ＊＊＊ xi, int T, int N)
{
int t;
for (t = 1; t <= T; t++)
free_dmatrix(xi[t], 1, N, 1, N);
xi ++;
free(xi);
}