hmm.cpp

经典的HMM算法的代码！以在文本中的应用为例

		for( int j=0 ; j<N ; j++ )
		{
			dMax = 0.0 ;
			dValue = 0.0 ;
			for( int i=0 ; i<N ; i++ )
			{
				dValue = delta[t-1][i] * A[i][j] ;
				if( dMax < dValue )
				{
					dMax = dValue ;
					Psi[t][j] = i ;
				}
			}
			delta[t][j] = dMax * B[j][O[t]] ;
		}
	}

	//step three : termination
	dMax = 0.0 ;
	int iIndex = 0 ;
	for( i=0 ; i<N ; i++ )
	{
		if( dMax < delta[T-1][i] )
		{
			dMax = delta[T-1][i] ;
			iIndex = i ;
		}
	}
	S[T-1] = iIndex ;
	for( t=T-2 ; t>=0 ; t-- )
		S[t] = Psi[t+1][S[t+1]] ;

	//free the memory of Psi
	iFreeMatrix( Psi, 0, T-1, 0, N-1 ) ;

	return dMax ;

}

double Hmm::Viterbi( int T, int* O, int*& S )
{
	//to check the observation sequence
	if( T == 0 )
		return FALSEVALUE ;
	if( delta != NULL )
	{
		//delta is not empty indicate that this is a new observation sequence
		//then reallocate the memory is needed
		//fisrt free the memory of beta
		dFreeMatrix( delta, 0, iPreT-1, 0, N-1 ) ;/*the fourth parament is no use here since the third one is always 0*/
		//set the iPreT with the current value
		iPreT = T ;
	}
	//allocate memory
	delta = dMatrix( 0, T-1, 0, N-1 ) ;

	//allocate the memory for the state sequence S
	S = (int*)malloc( T*sizeof(int) ) ;

	//***************** this is the viterbi algorithm********************************/

	//the matrix for back trace
	int** Psi ;
	//allocate memory for it
	Psi = iMatrix( 0, T-1, 0, N-1 ) ;

	//step one : initialization
	for( int i=0 ; i<N ; i++ )
	{
		delta[0][i] = pi[i] * B[i][O[0]] ;
		Psi[0][i] = 0 ;
	}

	//step two : induction
	double dMax = 0.0 ;
	double dValue = 0.0 ;
	for( int t=1 ; t<T ; t++ )
	{
		for( int j=0 ; j<N ; j++ )
		{
			dMax = 0.0 ;
			dValue = 0.0 ;
			for( int i=0 ; i<N ; i++ )
			{
				dValue = delta[t-1][i] * A[i][j] ;
				if( dMax < dValue )
				{
					dMax = dValue ;
					Psi[t][j] = i ;
				}
			}
			delta[t][j] = dMax * B[j][O[t]] ;
		}
	}

	//step three : termination
	dMax = 0.0 ;
	int iIndex = 0 ;
	for( i=0 ; i<N ; i++ )
	{
		if( dMax < delta[T-1][i] )
		{
			dMax = delta[T-1][i] ;
			iIndex = i ;
		}
	}
	S[T-1] = iIndex ;
	for( t=T-2 ; t>=0 ; t-- )
		S[t] = Psi[t+1][S[t+1]] ;

	//free the memory of Psi
	iFreeMatrix( Psi, 0, T-1, 0, N-1 ) ;

	return dMax ;

}

void Hmm::BaumWelch( int T, vector<int>& O, double& probInit, double& probFinal )
{
	//step one : calculate alpha and initial probability using forward procedure
	probInit = Forward( T, O ) ;

	//step two : calculate beta using backward procedure
	Backward( T, O ) ;//first it's equal to probInit

	//step three : allocate memory for gamma and calculate it
	gamma = ComputeGamma( T, N, probInit ) ;

	//step four : allocate memory for p[t][i][j] and calculate it
	//p[t][i][j] = p( Xt=i,Xt+1=j | O,u ), is the probability of transitting from state i to state j
	double*** p = NULL ;
	p = ComputeP( p, T, N, O, probInit ) ;

	//step five : recalculate all the parameters of hmm( pi, A and B, and also alpha,beta,gamma )
	//untill some consition is abtained
	RecomputeParameter( p, T, O, probInit, probFinal ) ;

	//step six : free all the memory that allocated before
	PFreeMatrix( p, T ) ;

}

void Hmm::BaumWelch( int T, int* O, double& probInit, double& probFinal )
{
	//step one : calculate alpha and initial probability using forward procedure
	probInit = Forward( T, O ) ;

	//step two : calculate beta using backward procedure
	Backward( T, O ) ;//first it's equal to probInit

	//step three : allocate memory for gamma and calculate it
	gamma = ComputeGamma( T, N, probInit ) ;

	//step four : allocate memory for p[t][i][j] and calculate it
	//p[t][i][j] = p( Xt=i,Xt+1=j | O,u ), is the probability of transitting from state i to state j
	double*** p = NULL ;
	p = ComputeP( p, T, N, O, probInit ) ;

	//step five : recalculate all the parameters of hmm( pi, A and B, and also alpha,beta,gamma )
	//untill some consition is abtained
	RecomputeParameter( p, T, O, probInit, probFinal ) ;

	//step six : free all the memory that allocated before
	PFreeMatrix( p, T ) ;
}

/********************************the core functions of Hmm end*******************************************/

/******************************** the function called by baumwelch ***************************************/
double** Hmm::ComputeGamma( int T, int N, double prob ) 
{
	//firstly allocate memory for gamma
	if( gamma )
	{
		//if gamma is not null, then we should free the memory that gamma ocuppied
		dFreeMatrix( gamma, 0, iPreT-1, 0, N-1 ) ;
	}
	gamma = dMatrix( 0, T-1, 0, N-1 ) ;

	//initialize
	for( int t=0 ; t<T ; t++ )
		for( int i=0 ; i<N ; i++ )
			gamma[t][i] = 0.0 ;

	//calculate gamma using alpha and beta
	double dDenominator = DENOMINATORLIMIT ;
	for( t=0 ; t<T ; t++ )
	{
		for( int i=0 ; i<N ; i++ )
		{
			gamma[t][i] = alpha[t][i] * beta[t][i] ;
			gamma[t][i] /= prob ;
		}
	}

	return gamma ;
}

double*** Hmm::ComputeP( double*** p ,int T, int N, vector<int>& O, double prob )
{
	//firstly allocate memory for p
	if( p )
	{
		//if p is not null, then we should free the memory that p ocuppied
		for( int t=0 ; t<T ; t++ )
			dFreeMatrix( p[t], 0, N-1, 0, N-1 ) ;
		delete[] p ;
	}
	p = new double** [T] ;
	for( int t=0 ; t<T ; t++ )
		p[t] = dMatrix( 0, N-1, 0, N-1 ) ;

	//initialize the 3-dimision matrix
	for( t=0 ; t<T ;t++ )
		for( int i=0 ; i<N ; i++ )
			for( int j=0 ; j<N ; j++ )
				p[t][i][j] = 0.0 ;

	//calculate it
	for( t=0 ; t<T-1 ; t++ ) 
	//attension : the time variable is from 0 to T-2, together t-1 time, because there are only t-1 transittions!
	{
		for( int i=0 ; i<N ; i++ )
		{
			for( int j=0 ; j<N ; j++ )
			{
				p[t][i][j] = alpha[t][i] * A[i][j] * B[j][O[t+1]] * beta[t+1][j] ;
				p[t][i][j] /= prob ;
			}
		}
	}

	return p ;
}

double*** Hmm::ComputeP( double*** p, int T, int N, int* O, double prob )
{
	//firstly allocate memory for p
	if( p )
	{
		//if p is not null, then we should free the memory that p ocuppied
		for( int t=0 ; t<T ; t++ )
			dFreeMatrix( p[t], 0, N-1, 0, N-1 ) ;
		delete[] p ;
	}
	p = new double** [T] ;
	for( int t=0 ; t<T ; t++ )
		p[t] = dMatrix( 0, N-1, 0, N-1 ) ;

	//initialize the 3-dimision matrix
	for( t=0 ; t<T ;t++ )
		for( int i=0 ; i<N ; i++ )
			for( int j=0 ; j<N ; j++ )
				p[t][i][j] = 0.0 ;

	//calculate it
	for( t=0 ; t<T-1 ; t++ ) 
	//attension : the time variable is from 0 to T-2, together t-1 time, because there are only t-1 transittions!
	{
		for( int i=0 ; i<N ; i++ )
		{
			for( int j=0 ; j<N ; j++ )
			{
				p[t][i][j] = alpha[t][i] * A[i][j] * B[j][O[t+1]] * beta[t+1][j] ;
				p[t][i][j] /= prob ;
			}
		}
	}

	return p ;
}

void Hmm::RecomputeParameter( double***& p, int T, vector<int>& O, double probInit, double& probFinal )
{
	//record the round number of do..while circle and jump out of it when certain number exceed
	int iRound = 0 ;

	double dRatio = 0.0 ;
	double probTemp = 0.0 ;
	double probOld = probInit ;
	do{
		//recalculate initial probability
		pi.clear() ;
		for( int i=0 ; i<N ; i++ )
			pi.push_back( gamma[0][i] ) ;

		//recalculate the transition matrix A probability
		//double dDenominatorA = 0.0 ;
		double dDenominatorA = DENOMINATORLIMIT ;
		double dNumeratorA = 0.0 ;
		for( i=0 ; i<N ; i++ )
		{
			//first calculate the denominator of A[i][j], that's the sum of gamma[t][i]
			//(excluded the time T-1 because there is no transition from the time T-1)
			dDenominatorA = DENOMINATORLIMIT ;
			for( int t=0 ; t<T-1 ; t++ )
				dDenominatorA += gamma[t][i] ;
			for( int j=0 ; j<N ; j++ )
			{
				//second calculate the numerator of A[i][j], that's the sum of p[t][i][j]
				//(excluded the time T-1 because there is no transition from the time T-1)
				dNumeratorA = 0.0 ;
				for( t=0 ; t<T-1 ; t++ )
					dNumeratorA += p[t][i][j] ;

				//recalculate the probability of A[i][j]
				A[i][j] = dNumeratorA / dDenominatorA ;
			}
		}

		//recalculate the emit probability matrix of B
		double dDenominatorB = DENOMINATORLIMIT ;
		double dNumeratorB = 0.0 ;
		for( int j=0 ; j<N ; j++ )
		{
			dDenominatorB = DENOMINATORLIMIT ;
			for( int t=0 ; t<T ; t++ )
				dDenominatorB += gamma[t][j] ;
			for( int k=0 ; k<M ; k++ )
			{
				dNumeratorB = 0.0 ;
				for( int t=0 ; t<T ;t++ )
				{
					if( O[t] == k )
						dNumeratorB += gamma[t][j] ;
				}

				//recalculate the probability of B[j][k]
				B[j][k] = dNumeratorB / dDenominatorB ;
			}
		}

		//calculate alpha and the probability once again, repeat the "step one"
		probTemp = Forward( T, O ) ;
		//calculate beta , repeat the "step two"
		Backward( T, O ) ;
		//allocate memory for gamma and calculate it, repeat the "step three"
		gamma = ComputeGamma( T, N, probTemp ) ;
		//recalculate p[t][i][j]
		p = ComputeP( p, T, N, O, probTemp ) ;

		dRatio = ( probTemp - probOld ) / probOld ;
		probOld = probTemp ;

		//when the circle number of do..while exceed ROUNDLIMIT times, break it
		//thus, we may not get stuck in a local maximum
		iRound++ ;
		if( iRound > ROUNDLIMIT )
			break ;

	}while( dRatio>RATIOLIMIT ) ;

	probFinal = probTemp ;

}

void Hmm::RecomputeParameter( double***& p, int T, int* O, double probInit, double& probFinal )
{
	//record the round number of do..while circle and jump out of it when certain number exceed
	int iRound = 0 ;

	double dRatio = 0.0 ;
	double probTemp = 0.0 ;
	double probOld = probInit ;
	do{
		//recalculate initial probability
		pi.clear() ;
		for( int i=0 ; i<N ; i++ )
			pi.push_back( gamma[0][i] ) ;

		//recalculate the transition matrix A probability
		//double dDenominatorA = 0.0 ;
		double dDenominatorA = DENOMINATORLIMIT ;
		double dNumeratorA = 0.0 ;
		for( i=0 ; i<N ; i++ )
		{
			//first calculate the denominator of A[i][j], that's the sum of gamma[t][i]
			//(excluded the time T-1 because there is no transition from the time T-1)
			dDenominatorA = DENOMINATORLIMIT ;
			for( int t=0 ; t<T-1 ; t++ )
				dDenominatorA += gamma[t][i] ;
			for( int j=0 ; j<N ; j++ )
			{
				//second calculate the numerator of A[i][j], that's the sum of p[t][i][j]
				//(excluded the time T-1 because there is no transition from the time T-1)
				dNumeratorA = 0.0 ;
				for( t=0 ; t<T-1 ; t++ )
					dNumeratorA += p[t][i][j] ;

				//recalculate the probability of A[i][j]
				A[i][j] = dNumeratorA / dDenominatorA ;
			}
		}

		//recalculate the emit probability matrix of B
		double dDenominatorB = DENOMINATORLIMIT ;
		double dNumeratorB = 0.0 ;
		for( int j=0 ; j<N ; j++ )
		{
			dDenominatorB = DENOMINATORLIMIT ;
			for( int t=0 ; t<T ; t++ )
				dDenominatorB += gamma[t][j] ;
			for( int k=0 ; k<M ; k++ )
			{
				dNumeratorB = 0.0 ;
				for( int t=0 ; t<T ;t++ )
				{
					if( O[t] == k )
						dNumeratorB += gamma[t][j] ;
				}

				//recalculate the probability of B[j][k]
				B[j][k] = dNumeratorB / dDenominatorB ;
			}
		}

		//calculate alpha and the probability once again, repeat the "step one"
		probTemp = Forward( T, O ) ;
		//calculate beta , repeat the "step two"
		Backward( T, O ) ;
		//allocate memory for gamma and calculate it, repeat the "step three"
		gamma = ComputeGamma( T, N, probTemp ) ;
		//recalculate p[t][i][j]
		p = ComputeP( p, T, N, O, probTemp ) ;

		dRatio = ( probTemp - probOld ) / probOld ;
		probOld = probTemp ;

		//when the circle number of do..while exceed ROUNDLIMIT times, break it
		//thus, we may not get stuck in a local maximum
		iRound++ ;
		if( iRound > ROUNDLIMIT )
			break ;

	}while( dRatio>RATIOLIMIT ) ;

	probFinal = probTemp ;
}

void Hmm::PFreeMatrix( double*** p, int T )
{
	for( int t=0 ; t<T ; t++ )
		dFreeMatrix( p[t], 0, N-1, 0, N-1 ) ;

	delete[] p ;
}
/****************** the function called by baumwelch end **********************************/

/************************the generating sequence function*******************************/
void Hmm::GenerateSequence( int iSeed, int T, vector<int>& O, vector<int>& S )
{
	//allocate memory for vector O and S
	O.resize( T ) ;
	S.resize( T ) ;

	//get the initial state
	S[0] = GetInitialState( iSeed ) ;
	O[0] = GetSymbol( S[0] ) ;

	//get the rest state and symbol
	for( int t=1 ; t<T ; t++ )
	{
		S[t] = GetNextState( S[t-1] ) ;
		O[t] = GetSymbol( S[t] ) ;
	}
}

void Hmm::GenerateSequence( int iSeed, int T, int* O, int*S )
{
	//get the initial state
	S[0] = GetInitialState( iSeed ) ;
	O[0] = GetSymbol( S[0] ) ;

	//get the rest state and symbol
	for( int t=1 ; t<T ; t++ )
	{
		S[t] = GetNextState( S[t-1] ) ;
		O[t] = GetSymbol( S[t] ) ;
	}
}

int Hmm::GetInitialState( int iSeed ) 
{
	double dVal = 0.0 ;
	double dAccum = 0.0 ;
	int iState = -1 ;

	//set the seed
	srand( iSeed ) ;
	//generate a value between 0 and 1
	dVal = rand() / RAND_MAX ;
	
	//get a state according the state probability
	for( int i=0 ; i<N ; i++ )
	{
		if( dVal < (pi[i]+dAccum) )
		{
			iState = i ;
			break ;
		}
		else
			dAccum += pi[i] ;
	}

	return iState ;
}

int Hmm::GetSymbol( int iState )
{
	double dVal = 0.0 ;
	double dAccum = 0.0 ;

	int iSymbol = -1 ;

	int iSeed = 0 ;
	//get the current time as the seed
	long ltime = 0 ;
	ltime = time( NULL ) ;
	iSeed = (unsigned)ltime/2 ;

	//set the seed
	srand( iSeed ) ;
	//generate a value between 0 and 1
	dVal = rand() / RAND_MAX ;

	//get an output symbol according the emition probability
	for( int k=0 ; k<M ; k++ )
	{
		if( dVal < (B[iState][k]+dAccum) )
		{
			iSymbol = k ;
			break ;
		}
		else
			dAccum += B[iState][k] ;
	}

	return iSymbol ;
}

int Hmm::GetNextState( int iState )
{
	double dVal = 0.0 ;
	double dAccum = 0.0 ;

	int iStateNext = -1 ;

	int iSeed = 0 ;
	//get the current time as the seed
	long ltime = 0 ;
	ltime = time( NULL ) ;
	iSeed = (unsigned)ltime/2 ;
	//set the seed
	srand( iSeed ) ;
	//generate a value between 0 and 1
	dVal = rand() / RAND_MAX ;

	//get the next state according the state transittion probability
	for( int i=0 ; i<M ; i++ )
	{
		if( dVal < (A[iState][i]+dAccum) )
		{
			iStateNext = i ;
			break ;
		}
		else
			dAccum += A[iState][i] ;
	}

	return iStateNext ;
}

/**********************the generating sequence function end **************************/