cheekdiscretehmm.cpp

读入图像序列

/***************************************************
*离散HMM的实现文件
*作者：周庚涛
*时间：2006.03.29
*第一次修改者：周庚涛
*第一次修改时间 ：2006.6.21
***************************************************/

#include "StdAfx.h"
#include "cheekdiscretehmm.h"
#include "math.h"

CCheekDiscreteHmm::CCheekDiscreteHmm()
{
	int i = 0;
	int j = 0;
	
	
	pi[0]=1.0;
	pi[1]=0.0;
	pi[2]=0.0;
	



	for(i=0; i<N; i++) //M个观察值
	{
		for(j=0; j<N; j++) 
		{
			if (i>j)
			{
				A[i][j]=0;
			}
		
		}
	}


	A[0][0]=A[1][1]= 10.0/11;
	A[2][2]=1.0;
	A[0][1]=1.0/22;
	A[0][2]=1.0/22;
	A[1][2]=1.0/11;

	
	for(i=0; i<N; i++) //M个观察值
	{
		for(j=0; j<M; j++) 
		{
			B[i][j] = 1.0/M;
		}
	}

}

CCheekDiscreteHmm::~CCheekDiscreteHmm()
{
}



void CCheekDiscreteHmm::ForwardWithScale(int O[T], double alpha[T][N],double scale[T], double & pprob)
{

	int i,j; 
	int t;

	double sum;
	/*
	 *	1 初始化
	 */
    scale[0] = 0.0;
	for(i=0; i<N; i++)
	{
		alpha[0][i] = pi[i] * B[i][O[0]];
		scale[0] += alpha[0][i];
	}
	for(i=0; i<N; i++)
	{
		alpha[0][i]/= scale[0];
	}

	/*
	 *	2递归
	 */

	for(t=0; t<T-1; t++)
	{
		scale[t+1] = 0.0;
		for(j=0; j<N; j++)
		{
			sum = 0.0;
			for(i=0; i<N;i++)
			{
				sum += alpha[t][i]*A[i][j];
			}
			alpha[t+1][j] = sum *B[j][O[t+1]];
			scale [t+1] += alpha[t+1][j];
		}
		for(j=0; j<N; j++)
		{
			alpha[t+1][j]/=scale[t+1];
		}
	}
	/*
	*	3 终止
	*/
	pprob = 0.0;
	for(t=0;t<T;t++)
	{
		pprob += log(scale[t]);

	}

	
}


/***************************************************************************
** 函数名称：BackwardWithScale
** 功能：后向算法估计参数（带比例因子修正）
** 参数：O：观察值序列
**       beta：运算中用到的临时数组
**       scale：比例因子数组
**       pprob：返回值，所要求的概率
**/
void CCheekDiscreteHmm::BackwardWithScale(int O[T], double beta[T][N],double scale[T]/*,double & pprob*/)
{
	int     i, j;   /* 状态指示 */
	int     t;      /* 时间下标 */
	double sum;
	
	
	/* 1. 初始化 */
	for (i = 0; i < N; i++)
	{
		beta[T-1][i] = 1.0/scale[T-1]; 
	}
	
	/* 2. 递归 */
	for (t = T - 2; t >= 0; t--) 
	{
		for (i = 0; i < N; i++) 
		{
			sum = 0.0;
			for (j = 0; j < N; j++)
			{
				sum += A[i][j] * (B[j][O[t+1]])*beta[t+1][j];
			}
			beta[t][i] = sum/scale[t];
		}
	}



}
/******************************************************************************
**函数名称：BaumWelch
**功能：BaumWelch算法
**参数：phmm：HMM模型指针
**      T：观察序列长度
**      O：观察序列
**      alpha，beta，gamma，pniter均为中间变量
**      plogprobinit：初始概率
**      最终概率
**/ 
int CCheekDiscreteHmm::BaumWelch(int O[NUM][T])
{
	int	i, j, k;
	int n;
	int	t, l = 0;
	int num=0;

	double alpha[T][N];
	double beta[T][N];
	double gamma[T][N];
	double xi[T-1][N][N];	
	double scale[T];

	double tempPro;
    double profirst=0.0, prosecond=0.0;
	
	double	logprobf;
	int state[T];
	double	numeratorA, denominatorA;
	double	numeratorB, denominatorB;	
	
	double  deltaprev;
	
	deltaprev = 10e-70;	
	
	do  
	{	
		for (n=0; n< NUM; n++)
		{
			ViterbiLog(O[n],state,tempPro);
			profirst += tempPro;			
		}
	

		
		/* 重新估计转移矩阵*/
		for (i = 0; i < N; i++) 
		{ 
			denominatorA = 0.0;
			for(n=0; n<NUM; n++)
			{
				ForwardWithScale(O[n], alpha, scale, logprobf);
				BackwardWithScale(O[n], beta, scale/*, logprobb*/);
				ComputeGamma(alpha, beta, gamma);
				ComputeXi(O[n], alpha, beta, xi);
				for (t = 0; t <  T - 1; t++) 
				{
					denominatorA += gamma[t][i];
				}
			}
						
			for (j = i; j < N; j++) 
			{
				numeratorA = 0.0;
				for(n=0; n<NUM; n++)
				{
					ForwardWithScale(O[n], alpha, scale, logprobf);
					BackwardWithScale(O[n], beta, scale/*, logprobb*/);
					ComputeGamma(alpha, beta, gamma);
					ComputeXi(O[n], alpha, beta, xi);
					for (t = 0; t <  T - 1; t++) 
					{
						numeratorA += xi[t][i][j];
					}					
				}		
			
				A[i][j] = 0.001 + 0.999*numeratorA/denominatorA;
			}
			//?/////////////////////////////////////

            /*观察矩阵 */
			denominatorB = 0.0;
			for(n=0; n<NUM; n++)
			{
				ForwardWithScale(O[n], alpha, scale, logprobf);
				BackwardWithScale(O[n], beta, scale/*, logprobb*/);
				ComputeGamma(alpha, beta, gamma);
				ComputeXi(O[n], alpha, beta, xi);
				for (t = 0; t <  T ; t++) 
				{
					denominatorB += gamma[t][i];
				}
			}
			//denominatorB = denominatorA + gamma[T-1][i]; 
			for (k = 0; k < M; k++) 
			{
				numeratorB = 0.0;
				for(n=0; n<NUM; n++)
				{
					ForwardWithScale(O[n], alpha, scale, logprobf);
					BackwardWithScale(O[n], beta, scale/*, logprobb*/);
					ComputeGamma(alpha, beta, gamma);
					ComputeXi(O[n], alpha, beta, xi);
				
					for (t = 0; t < T; t++) 
					{
						if (O[n][t] == k) 
						{
							numeratorB += gamma[t][i];
						}
					}
				}
				
				B[i][k] = 0.001 + 0.999*numeratorB/denominatorB;
			}
		}
		for (n=0; n< NUM; n++)
		{
			ViterbiLog(O[n],state,tempPro);
			prosecond += tempPro;			
		}		
		num++;	
 	}while (fabs((prosecond-profirst)/prosecond)>exp(-4));
	/* 如果差的不太大，表明收敛，退出 */ 
	return 0;
	
	

}


void CCheekDiscreteHmm::ComputeGamma(double alpha[T][N], double beta[T][N], double gamma[T][N])
{
	int i, j;
	int	t;
	double	denominator;
	
	for (t = 0; t < T; t++) 
	{
		denominator = 0.0;
		for (j = 0; j < N; j++) 
		{
			gamma[t][j] = alpha[t][j]*beta[t][j];
			denominator += gamma[t][j];
		}
		
		for (i = 0; i < N; i++) 
		{
			gamma[t][i] = gamma[t][i]/denominator;
		}
	}
}
void CCheekDiscreteHmm::ComputeXi(int O[T], double alpha[T][N], double beta[T][N], double xi[T][N][N])
{
	int i, j;
	int t;
	double sum;
	
	for (t = 0; t < T - 1; t++) 
	{
		sum = 0.0;	
		for (i = 0; i < N; i++) 
		{
			for (j = 0; j < N; j++) 
			{
				xi[t][i][j] = alpha[t][i]*beta[t+1][j]* A[i][j] * B[j][O[t+1]];
				sum += xi[t][i][j];
			}
		}
		for (i = 0; i < N; i++) 
		{
			for (j = 0; j < N; j++)
			{
				xi[t][i][j]  /= sum;
			}
		}
	}

}



/**************************************************************************
** 函数名称：ViterbiLog
** 功能：Viterbi算法
** 参数：phmm：HMM结构指针
**       T：观察值的个数
**       O：观察序列
**       delta，psi为中间变量
**       q：求得的最佳状态序列
**       pprob：概率
*****************************************************************************/
void CCheekDiscreteHmm::ViterbiLog(int O[T], int q[T], double & pprob)
{
	int     i, j;   /* 状态下标 */
	int     t;      /* 时间下标 */

	double delta[T][N]; 
	int psi[T][N];
	
	int     maxvalind;
	double  maxval, val;
	double  biot[N][T];
	double  PITemp[N];
	double ATemp[N][N];
	double inf = 999999;
	
	/* 0. 预处理 */
	
	for (i = 0; i < N; i++) 
	{
		if (pi[i] != 0)
		{
			PITemp[i] = log(pi[i]);
		}
		else
		{
			PITemp[i] = -inf;

		}
	}
	for (i = 0; i < N; i++) 
	{
		for (j = 0; j < N; j++) 
		{
			if (A[i][j] != 0)
			{
				ATemp[i][j] = log(A[i][j]);

			}
			else
			{
				ATemp[i][j] = -inf;

			}
		}
	}
	
	//	biot = dmatrix(1, N, 1, T);////////////////////////////
	for (i = 0; i < N; i++) 
	{
		for (t = 0; t < T; t++) 
		{
			if (B[i][O[t]] != 0) 
			{
				biot[i][t] = log(B[i][O[t]]);
			}
			else
			{
				biot[i][t] = -inf;

			}
		}
	}
	
	/* 1. 初始化  */
	
	for (i = 0; i < N; i++) 
	{
		delta[0][i] =PITemp[i] + biot[i][0];
		psi[0][i] = 0;
	}
	
	/* 2. 递归 */
	
	for (t = 1; t < T; t++) 
	{
		for (j = 0; j < N; j++) 
		{
			maxval = -VITHUGE;
			maxvalind = 1;
			for (i = 0; i < N; i++) 
			{
				val = delta[t-1][i] + (ATemp[i][j]);
				if (val > maxval) 
				{
					maxval = val;
					maxvalind = i;
				}
			}
			
			delta[t][j] = maxval + biot[j][t]; 
			psi[t][j] = maxvalind;
			
		}
	}
	
	/* 3. 终止 */
	
	pprob = -VITHUGE;
	q[T-1] = 1;
	for (i = 0; i < N; i++) 
	{

		if (delta[T-1][i] > pprob) 
		{
			pprob = delta[T-1][i];
			q[T-1] = i;
		}
	}

	
	/* 4. 回溯 */
	
	for (t = T - 2; t >= 0; t--)
	{
		q[t] = psi[t+1][q[t+1]];	
	}
}