<P>/*<BR>** File:Viterbi.cpp<BR>** 功能:给定HMM和观察序列,求最可能的状态<BR>*/</P>
<P>//#include "StdAfx.h"<BR>#include <math.h><BR>#include "hmm.h"<BR>#include "nrutil.h"</P>
<P>#define VITHUGE 100000000000.0</P>
<P>/**************************************************************************<BR>** 函数名称:Viterbi<BR>** 功能:Viterbi算法<BR>** 参数:phmm:HMM结构指针<BR>** T:观察值的个数<BR>** O:观察序列<BR>** delta,psi为中间变量<BR>** q:求得的最佳状态序列<BR>** pprob:概率<BR>**/<BR>void Viterbi(HMM *phmm, int T, int *O, double **delta, int **psi, <BR> int *q, double *pprob)<BR>{<BR> int i, j; /* 状态下标 */<BR> int t; /* 时间下标 */ </P>
<P> int maxvalind;<BR> double maxval, val;</P>
<P> /* 1. 初始化 */<BR> <BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> delta[1] = phmm->pi * (phmm->B[O[1]]);<BR> psi[1] = 0;<BR> } </P>
<P> /* 2. 递归 */<BR> <BR> for (t = 2; t <= T; t++) <BR> {<BR> for (j = 1; j <= phmm->N; j++) <BR> {<BR> maxval = 0.0;<BR> maxvalind = 1; <BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> val = delta[t-1]*(phmm->A[j]);<BR> if (val > maxval) {<BR> maxval = val; <BR> maxvalind = i; <BR> }<BR> }<BR> <BR> delta[t][j] = maxval*(phmm->B[j][O[t]]);<BR> psi[t][j] = maxvalind; </P>
<P> }<BR> }</P>
<P> /* 3. 终止 */</P>
<P> *pprob = 0.0;<BR> q[T] = 1;<BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> if (delta[T] > *pprob) <BR> {<BR> *pprob = delta[T]; <BR> q[T] = i;<BR> }<BR> }</P>
<P> /* 4. Path (state sequence) backtracking */</P>
<P> for (t = T - 1; t >= 1; t--)<BR> q[t] = psi[t+1][q[t+1]];</P>
<P>}</P>
<P>/**************************************************************************<BR>** 函数名称:ViterbiLog<BR>** 功能:Viterbi算法<BR>** 参数:phmm:HMM结构指针<BR>** T:观察值的个数<BR>** O:观察序列<BR>** delta,psi为中间变量<BR>** q:求得的最佳状态序列<BR>** pprob:概率<BR>**/<BR>void ViterbiLog(HMM *phmm, int T, int *O, double **delta, int **psi,<BR> int *q, double *pprob)<BR>{<BR> int i, j; /* 状态下标 */<BR> int t; /* 时间下标 */<BR> <BR> int maxvalind;<BR> double maxval, val;<BR> double **biot;</P>
<P> /* 0. 预处理 */</P>
<P> for (i = 1; i <= phmm->N; i++) <BR> phmm->pi = log(phmm->pi);<BR> for (i = 1; i <= phmm->N; i++) <BR> for (j = 1; j <= phmm->N; j++) <BR> {<BR> phmm->A[j] = log(phmm->A[j]);<BR> }</P>
<P> biot = dmatrix(1, phmm->N, 1, T);<BR> for (i = 1; i <= phmm->N; i++) <BR> for (t = 1; t <= T; t++) <BR> {<BR> biot[t] = log(phmm->B[O[t]]);<BR> }<BR> <BR> /* 1. 初始化 */<BR> <BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> delta[1] = phmm->pi + biot[1];<BR> psi[1] = 0;<BR> }<BR> <BR> /* 2. 递归 */<BR> <BR> for (t = 2; t <= T; t++) <BR> {<BR> for (j = 1; j <= phmm->N; j++) <BR> {<BR> maxval = -VITHUGE;<BR> maxvalind = 1;<BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> val = delta[t-1] + (phmm->A[j]);<BR> if (val > maxval) <BR> {<BR> maxval = val;<BR> maxvalind = i;<BR> }<BR> }<BR> <BR> delta[t][j] = maxval + biot[j][t]; <BR> psi[t][j] = maxvalind;<BR> <BR> }<BR> }<BR> <BR> /* 3. 终止 */<BR> <BR> *pprob = -VITHUGE;<BR> q[T] = 1;<BR> for (i = 1; i <= phmm->N; i++) <BR> {<BR> if (delta[T] > *pprob) <BR> {<BR> *pprob = delta[T];<BR> q[T] = i;<BR> }<BR> }<BR> <BR> <BR> /* 4. 回溯 */</P>
<P> for (t = T - 1; t >= 1; t--)<BR> q[t] = psi[t+1][q[t+1]];</P>
<P>}<BR></P> |