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);
}