typedef struct
{
int N; /* 隐藏状态数目;Q={1,2,…,N} */
int M; /* 观察符号数目; V={1,2,…,M}*/
double **A; /* 状态转移矩阵A[1..N][1..N]. a[i][j] 是从t时刻状态i到t+1时刻状态j的转移概率 */
double **B; /* 混淆矩阵B[1..N][1..M]. b[j][k]在状态j时观察到符合k的概率。*/
double *pi; /* 初始向量pi[1..N],pi[i] 是初始状态概率分布 */
} HMM;前向算法程序示例如下:
/*函数参数说明:*phmm:已知的HMM模型;T:观察符号序列长度;*O:观察序列;**alpha:局部概率;*pprob:最终的观察概率
*/
void Forward(HMM *phmm, int T, int *O, double **alpha, double *pprob)
{int i, j; /* 状态索引 */int t; /* 时间索引 */double sum; /*求局部概率时的中间值 *//* 1. 初始化:计算t=1时刻所有状态的局部概率: */for (i = 1; i <= phmm->N; i++)alpha[1][i] = phmm->pi[i]* phmm->B[i][O[1]];/* 2. 归纳:递归计算每个时间点,t=2,… ,T时的局部概率 */for (t = 1; t < T; t++){for (j = 1; j <= phmm->N; j++){sum = 0.0;for (i = 1; i <= phmm->N; i++)sum += alpha[t][i]* (phmm->A[i][j]);alpha[t+1][j] = sum*(phmm->B[j][O[t+1]]);}}/* 3. 终止:观察序列的概率等于T时刻所有局部概率之和*/*pprob = 0.0;for (i = 1; i <= phmm->N; i++)*pprob += alpha[T][i];
}维比特算法:
void Viterbi(HMM *phmm, int T, int *O, double **delta, int **psi,int *q, double *pprob)
{int i, j; /* state indices */int t; /* time index */int maxvalind;double maxval, val;/* 1. Initialization */for (i = 1; i <= phmm->N; i++){delta[1][i] = phmm->pi[i] * (phmm->B[i][O[1]]);psi[1][i] = 0;}/* 2. Recursion */for (t = 2; t <= T; t++){for (j = 1; j <= phmm->N; j++){maxval = 0.0;maxvalind = 1;for (i = 1; i <= phmm->N; i++){val = delta[t-1][i]*(phmm->A[i][j]);if (val > maxval){maxval = val;maxvalind = i;}}delta[t][j] = maxval*(phmm->B[j][O[t]]);psi[t][j] = maxvalind;}}/* 3. Termination */*pprob = 0.0;q[T] = 1;for (i = 1; i <= phmm->N; i++){if (delta[T][i] > *pprob){*pprob = delta[T][i];q[T] = i;}}/* 4. Path (state sequence) backtracking */for (t = T – 1; t >= 1; t–)q[t] = psi[t+1][q[t+1]];
}后向算法:
void Backward(HMM *phmm, int T, int *O, double **beta, double *pprob)
{int i, j; /* state indices */int t; /* time index */double sum;/* 1. Initialization */for (i = 1; i <= phmm->N; i++)beta[T][i] = 1.0;/* 2. Induction */for (t = T - 1; t >= 1; t--){for (i = 1; i <= phmm->N; i++){sum = 0.0;for (j = 1; j <= phmm->N; j++)sum += phmm->A[i][j] *(phmm->B[j][O[t+1]])*beta[t+1][j];beta[t][i] = sum;}}/* 3. Termination */*pprob = 0.0;for (i = 1; i <= phmm->N; i++)*pprob += beta[1][i];
}前向后向算法:
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, threshold;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 |intial model) */BackwardWithScale(phmm, T, O, beta, scale, &logprobb);ComputeGamma(phmm, T, alpha, beta, gamma);ComputeXi(phmm, T, O, alpha, beta, xi);logprobprev = logprobf;do{/* reestimate frequency of state i in time t=1 */for (i = 1; i <= phmm->N; i++)phmm->pi[i] = .001 + .999*gamma[1][i];/* reestimate transition matrix and symbol prob ineach state */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);/* compute difference between log probability oftwo iterations */delta = logprobf - logprobprev;logprobprev = logprobf;l++;}while (delta > DELTA); /* if log probability does notchange much, exit */*pniter = l;*plogprobfinal = logprobf; /* log P(O|estimated model) */FreeXi(xi, T, phmm->N);free_dvector(scale, 1, T);
}