# HMM前向算法，维比特算法，后向算法，前向后向算法代码

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 in
each 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 of
two iterations */
delta = logprobf - logprobprev;
logprobprev = logprobf;
l++;
}
while (delta > DELTA); /* if log probability does not
change much, exit */
*pniter = l;
*plogprobfinal = logprobf; /* log P(O|estimated model) */
FreeXi(xi, T, phmm->N);
free_dvector(scale, 1, T);
}