/***      File:   baumwelch.cpp**      功能：根據給定的觀察序列，用BaumWelch算法估計HMM模型參數*/
#include "StdAfx.h"#include <stdio.h> #include "nrutil.h"#include "hmm.h"#include <math.h>#define DELTA 0.001

/******************************************************************************
**函數名稱：BaumWelch
**功能：BaumWelch算法
**參數：phmm：HMM模型指針
**      T：觀察序列長度
**      O：觀察序列
**      alpha，beta，gamma，pniter均為中間變量
**      plogprobinit：初始概率
**      最終概率
**/ 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; /* log P(O |初始狀態) */	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);}