hmmutils.c

matlab 例程

	/*状态序列返回*/
	for(t=T-1;t>=1;t--)
		q[T]=psi[t+1][q[t+1]];

}
/*********************************
函数：ViterbiLog.c
功能：给定HMM和观察值序列，求最可能的状态
参数：phmm:HMM结构指针
      T:观察值的个数
	  O：观察值序列
	  deta,psi:中间变量
	  q:求的最佳状态序列
	  pprob:概率
*********************************/


void ViterbiLog(struct HMM *phmm,int T,int **O,double **delta,int **psi,int **q,double *pprob,int k)
{
	int i,j;
	int t;
	int maxvalind;
	double maxval,val;
	double **biot;
    struct HMM *phmm1;
   phmm1=(struct HMM*)malloc(sizeof(struct HMM*));
phmm1->pi=(double *)dvector(1,phmm->N);
    phmm1->A=(double **)dmatrix(1,phmm->N,1,phmm->N);
phmm1->B=(double **)dmatrix(1,phmm->N,1,phmm->M);
	/*预处理*/
	for(i=1;i<=phmm->N;i++)
		phmm1->pi[i]=log(phmm->pi[i]);
   	for(i=1;i<=phmm->N;i++)
       	for(j=1;j<=phmm->N;j++)
        {
			phmm1->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][k]]);
		 }
	/*初始化*/
	for(i=1;i<=phmm->N;i++)
	{
		delta[1][i]=phmm1->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]+(phmm1->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][k]=1;
	for(i=1;i<=phmm->N;i++)
	{
		if(delta[T][i]>*pprob)
		{
			*pprob=delta[T][i];
			q[T][k]=i;
		}
	}
	/*状态序列返回*/
	for(t=T-1;t>=1;t--)
		q[t][k]=psi[t+1][q[t+1][k]];

}
/***********************
函数BaumWelch
功能：BaumWelch算法
参数：phmm:HMM模型指针
      T:观察序列的长度
	  O：观察序列
	  alpha,beta,gamma,pniter均为中间变量
	  plogprobinit:初始概率
	  plogprobfinal:最终概率
*****************************/
#define DELTA -10e-1
void BaumWelch(struct HMM *phmm,int T,int **O,double **alpha,double **beta,
			   double **gamma,int *pniter,double *plogprobinit,double *plogprobfinal,int x,int y)
{
	int i,j;
	int l=1,t;
	int k;
	double logprobf,logprobb;
	double *numeratorA,*denominatorA;

	double sum1=0.0;
	double sum2=0.0;
    double sum3=0.0;
	double sum4=0.0;
    double sum5=0.0;
	double *numeratoru,*numerators,*denominatorB;
	double ***xi,*scale;
	double delta,deltaprev,logprobprev;
	deltaprev=10e-7;
	xi=AllocXi(T,2);
	scale=dvector(l,T);

numeratorA=dvector(l,(x/8)*(y/8));
denominatorA=dvector(l,(x/8)*(y/8));
numeratoru=dvector(l,(x/8)*(y/8));
numerators=dvector(l,(x/8)*(y/8));
denominatorB=dvector(l,(x/8)*(y/8));
logprobprev=logprobf=10;
do{	
for(i=1;i<=phmm->N;i++)
{ for(j=1;j<=phmm->N;j++)
 {

	for(k=1;k<=(x/8)*(y/8);k++)
	{
    
	ForwardWithScale(phmm,T,O,alpha,scale,&logprobf,k);
	*plogprobinit=logprobf;/*logP(O|初始状态)*/
	BackwardWithScale(phmm,T,O,beta,scale,&logprobb,k);
	ComputeGamma(phmm,T,alpha,beta,gamma);
	ComputeXi(phmm,T,O,alpha,beta,xi,k);
	
    
     /*重新估计t=1时，状态i的概率*/
		
			phmm->pi[i]+=.001+.009*gamma[1][i];
		    
		   
		
	/*从新估计转移矩阵和观察矩阵*/
	
        
         denominatorA[k]=0.0;
		 for(t=1;t<T;t++)
			 denominatorA[k]+=gamma[t][i];
		 
		 
			 numeratorA[k]=0.0;
			for(t=1;t<T;t++)
               numeratorA[k]+=xi[t][i][j];
			
		 
         
		
	
      denominatorB[k]=denominatorA[k]+gamma[T][i];
		 numeratoru[k]=0.0;
		for(t=1;t<=T;t++)
			 numeratoru[k]+=gamma[t][i]*O[t][k];
			 numerators[k]=0.0; 
       for(t=1;t<=T;t++)
			  numerators[k]+=gamma[t][i]*(O[t][k]-phmm->u[i])*(O[t][k]-phmm->u[i]);
			 
}
 	
    phmm->pi[i]/=((x/8)*(y/8));
    

	    for(k=1;k<=(x/8)*(y/8);k++)
	{            sum1+=numeratorA[k];
	             sum2+=denominatorA[k];
                 sum3+=numeratoru[k];
	             sum4+=denominatorB[k];
				 sum5+=numerators[k];
	}
			     phmm->A[i][j]=.001+.999*sum1/sum2;
                 phmm->u[i]=.001+.999*sum3/sum4;
	             phmm->sigma[i]=.001+.999*sum5/sum4;
		
}
}         
		 
         
		

	/*计算两次的概率差*/
	delta=logprobf-logprobprev;
		logprobprev=logprobf;
}	
	
	while(delta<DELTA);
	*pniter=1;
	*plogprobfinal=logprobf;/*logP(O|estimated model)*/
    FreeXi(xi,T,phmm->N);
	free_dvector(scale,l,T);
}

void BaumWelch0(struct HMM *phmm,int T,int **O,double **alpha,double **beta,
			   double **gamma,int *pniter,double *plogprobinit,double *plogprobfinal,int x,int y)
{
	int i,j;
	int l=1,t;
	int k;
	double logprobf,logprobb;
	double *numeratorA,*denominatorA;

	double sum1=0.0;
	double sum2=0.0;
    double sum3=0.0;
	double sum4=0.0;
    double sum5=0.0;
	double *numeratoru,*numerators,*denominatorB;
	double ***xi,*scale;
	double delta,deltaprev,logprobprev;
	deltaprev=10e-7;
	xi=AllocXi(T,2);
	scale=dvector(l,T);

numeratorA=dvector(l,x*y);
denominatorA=dvector(l,x*y);
numeratoru=dvector(l,x*y);
numerators=dvector(l,x*y);
denominatorB=dvector(l,x*y);
logprobprev=logprobf=10;
do{	
for(i=1;i<=phmm->N;i++)
{ for(j=1;j<=phmm->N;j++)
 {

	for(k=1;k<=x*y;k++)
	{
    
	ForwardWithScale(phmm,T,O,alpha,scale,&logprobf,k);
	*plogprobinit=logprobf;/*logP(O|初始状态)*/
	BackwardWithScale(phmm,T,O,beta,scale,&logprobb,k);
	ComputeGamma(phmm,T,alpha,beta,gamma);
	ComputeXi(phmm,T,O,alpha,beta,xi,k);
	
    
     /*重新估计t=1时，状态i的概率*/
		
			phmm->pi[i]+=.001+.009*gamma[1][i];
		    
		   
		
	/*从新估计转移矩阵和观察矩阵*/
	
        
         denominatorA[k]=0.0;
		 for(t=1;t<T;t++)
			 denominatorA[k]+=gamma[t][i];
		 
		 
			 numeratorA[k]=0.0;
			for(t=1;t<T;t++)
               numeratorA[k]+=xi[t][i][j];
			
		 
         
		
	
      denominatorB[k]=denominatorA[k]+gamma[T][i];
		 numeratoru[k]=0.0;
		for(t=1;t<=T;t++)
			 numeratoru[k]+=gamma[t][i]*O[t][k];
			 numerators[k]=0.0; 
       for(t=1;t<=T;t++)
			  numerators[k]+=gamma[t][i]*(O[t][k]-phmm->u[i])*(O[t][k]-phmm->u[i]);
			 
}
 	
    phmm->pi[i]/=(x*y);
    

	    for(k=1;k<=x*y;k++)
	{            sum1+=numeratorA[k];
	             sum2+=denominatorA[k];
                 sum3+=numeratoru[k];
	             sum4+=denominatorB[k];
				 sum5+=numerators[k];
	}
			     phmm->A[i][j]=.001+.999*sum1/sum2;
                 phmm->u[i]=.001+.999*sum3/sum4;
	             phmm->sigma[i]=.001+.999*sum5/sum4;
		
}
}         
		 
         
		

	/*计算两次的概率差*/
	delta=logprobf-logprobprev;
		logprobprev=logprobf;
}	
	
	while(delta<DELTA);
	*pniter=1;
	*plogprobfinal=logprobf;/*logP(O|estimated model)*/
    FreeXi(xi,T,phmm->N);
	free_dvector(scale,l,T);
}

/*非用户接口函数*/
void ComputeGamma(struct 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(struct HMM *phmm,int T,int **O,double **alpha,double **beta,
			   double ***xi,int k)
{
    int i,j;
	int t;
	double sum;
	for(t=1;t<T;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][k]]);
			   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);
}