matlab中hmm的使用示例

刚开始用matlab中的hmm函数的童鞋想必都对其参数或者返回值或多或少会有一些疑问，matlab官方文档 http://www.mathworks.cn/cn/help/stats/hidden-markov-models-hmm.html 对此的解释又是相当简洁，本文写了一些博主对这些函数的理解，如果有什么疑问或者指正的话欢迎提出O(∩_∩)O~

1. hmmdecode

1.1 matlab的运行

>> t = [.7, .3; .4, .6]
t =   %transmission matrix
    0.7000    0.3000
    0.4000    0.6000

>> e = [.1, .1 .8; .8, .1, .1]
e =  %emission matrix
    0.1000    0.1000    0.8000
    0.8000    0.1000    0.1000

>> s = [1, 2, 3]  %sequence
s =
     1     2     3

>> [p,logPseq,fs,bs,scale] = hmmdecode(s,t,e)
p =  %calculates the posterior state probabilities of the sequence seq
    0.2488    0.5771    0.9039
    0.7512    0.4229    0.0961

logPseq =  %the logarithm of the probability of sequence seq
   -4.2114

fs =  %fs and bs returns forward and backward probabilities of the sequence scaled by S.
    1.0000    0.2258    0.4677    0.9039
         0    0.7742    0.5323    0.0961

bs =
    1.0000    1.1020    1.2337    1.0000
    1.6445    0.9703    0.7946    1.0000

scale =
    1.0000    0.3100    0.1000    0.4782

1.2 返回值验证

写程序对返回值进行验证，程序见文后附录。

首先看程序运行后前向法和后向法的输出值，此处的fwd和bwd分别为前、后向法的先验概率。

验算可知matlab中的fs，bs为先验概率经过scale后的结果：

fwd：0.9039*0.31*0.1*0.4782 = 0.0134

bwd：1.102*0.4782*0.1 = 0.0527

此外，可知logPseq为sequence出现概率的log值，为什么返回的是log值还请大神赐教了~

>> exp(-4.2114)

ans =

    0.0148

并且，fwd[0,2]+fwd[1,2]=0.0148

为什么decode参数中木有初始状态矩阵？初始状态矩阵在matlab中被处理为[1,0]与transmission matix的乘积。

2.hmmtrain

>> [guessTR,guessE,logliks] = hmmtrain(s,t,e)

guessTR =
    0.0000    1.0000
    1.0000    0.0000

guessE =
    0.0000    1.0000    0.0000
    0.5000    0.0000    0.5000

logliks =
  Columns 1 through 9 
   -4.2114   -3.2519   -3.1930   -3.1534   -3.0951   -2.9867   -2.7676   -2.3981   -1.9886

  Columns 10 through 14 
   -1.6145   -1.4153   -1.3867   -1.3863   -1.3863

训练14次后收敛，结束后seq在此hmm中出现的概率为：

>> exp(-1.3863)
ans =
    0.2500

结果与代码运行相同，可能是由于matlab中初始状态矩阵做了特殊处理的原因，收敛的循环次数稍稍有点区别，但程序运行速度上matlab还是快了非常多~

附录：

hmm 代码(C#)， 其中train用了baum-welch算法

public void train(IEnumerable<int[]> trainsequence)    {        double[] pi_new = new double[numStates];        double[,] a_new = new double[numStates, numStates];        double[,] b_new = new double[numStates, numObservations];        double pi_count = 0.0f;        for (int i = 0; i < numStates; i++)        {            double numerator = 0;            double denominator = 0;            foreach (int[] sequence in trainsequence)            {                double[,] fwd = this.forwardProc(sequence);                double[,] bwd = this.backwardProc(sequence);                double[,] gamma_numerator = new double[sequence.Length, numStates];                double[] gamma_denominator = new double[sequence.Length];                double prob = this.getProbability(sequence);                for (int ii = 0; ii < numStates; ii++)                {                    gamma_numerator[0, ii] = fwd[ii, 0] * bwd[ii, 0];                    gamma_denominator[0] += gamma_numerator[0, ii];                }                numerator += gamma_numerator[0, i];                denominator += gamma_denominator[0];            }            pi_new[i] = numerator / denominator;            pi_count += pi_new[i];        }        for (int i = 0; i < numStates; i++)        {            for (int j = 0; j < numStates; j++)            {                double numerator = 0;                double denominator = 0;                foreach (int[] sequence in trainsequence)                {                    double[,] fwd = this.forwardProc(sequence);                    double[,] bwd = this.backwardProc(sequence);                    double[, ,] xi_numerator = new double[sequence.Length - 1, numStates, numStates];                    double[,] gamma_numerator = new double[sequence.Length, numStates];                    double[] xi_denominator = new double[sequence.Length - 1];                    double[] gamma_denominator = new double[sequence.Length];                    double prob = this.getProbability(sequence);                    double numerator_innersum = 0;                    double denominator_innersum = 0;                    for (int t = 0; t <= sequence.Length - 1; t++)                    {                        for (int ii = 0; ii < numStates; ii++)                        {                            gamma_numerator[t, ii] = fwd[ii, t] * bwd[ii, t];                            gamma_denominator[t] += gamma_numerator[t, ii];                            if (t != sequence.Length - 1)                            {                                for (int jj = 0; jj < numStates; jj++)                                {                                    xi_numerator[t, ii, jj] = fwd[ii, t] * TransitionProbabilities[ii, jj] * EmissionProbabilities[jj, sequence[t + 1]] * bwd[jj, t + 1];                                    xi_denominator[t] += xi_numerator[t, ii, jj];                                }                            }                        }                    }                    for (int t = 0; t < sequence.Length - 1; t++)                    {                        numerator_innersum += xi_numerator[t, i, j] / xi_denominator[t];                        denominator_innersum += gamma_numerator[t, i] / gamma_denominator[t];                    }                    numerator += numerator_innersum / prob;                    denominator += denominator_innersum / prob;                }                a_new[i, j] = numerator / denominator;            }        }        // recalculate emission probability b        for (int i = 0; i < numStates; i++)        {            for (int j = 0; j < numObservations; j++)            {                double numerator = 0;                double denominator = 0;                foreach (int[] sequence in trainsequence)                {                    double[,] fwd = this.forwardProc(sequence);                    double[,] bwd = this.backwardProc(sequence);                    double prob = this.getProbability(sequence);                    double[,] gamma_numerator = new double[sequence.Length, numStates];                    double[] gamma_denominator = new double[sequence.Length];                    double numerator_innersum = 0;                    double denominator_innersum = 0;                    for (int t = 0; t <= sequence.Length - 1; t++)                    {                        for (int ii = 0; ii < numStates; ii++)                        {                            gamma_numerator[t, ii] = fwd[ii, t] * bwd[ii, t];                            gamma_denominator[t] += gamma_numerator[t, ii];                        }                    }                    for (int t = 0; t <= sequence.Length - 1; t++)                    {                        if (sequence[t] == j)                        {                            numerator_innersum += gamma_numerator[t, i] / gamma_denominator[t];                        }                        denominator_innersum += gamma_numerator[t, i] / gamma_denominator[t];                    }                    numerator += (1 / prob) * numerator_innersum;                    denominator += (1 / prob) * denominator_innersum;                }                b_new[i, j] = numerator / denominator;            }        }        this.InitialStateProbabilities = pi_new;        this.TransitionProbabilities = a_new;        this.EmissionProbabilities = b_new;    }    protected double[,] forwardProc(int[] o)    {        double[,] f = new double[numStates, o.Length];        for (int l = 0; l < numStates; l++)        {            f[l, 0] = InitialStateProbabilities[l] * EmissionProbabilities[l, o[0]];            //Debug.Log (string.Format ("f{0}, 0 : {1} {2} {3}", l, f[l, 0], InitialStateProbabilities[l], EmissionProbabilities[l, o[0]]));        }        for (int i = 1; i < o.Length; i++)        {            for (int k = 0; k < numStates; k++)            {                double sum = 0;                for (int l = 0; l < numStates; l++)                {                    sum += f[l, i - 1] * TransitionProbabilities[l, k];                }                f[k, i] = sum * EmissionProbabilities[k, o[i]];                //Debug.Log (string.Format ("f{0}, {1} : {2} ", k, i, f[k, i]));            }        }        return f;    }        protected double[,] backwardProc(int[] o)    {        int T = o.Length;        double[,] bwd = new double[numStates, T];        /* Basisfall */        for (int i = 0; i < numStates; i++)            bwd[i, T - 1] = 1;        /* Induktion */        for (int t = T - 2; t >= 0; t--)        {            for (int i = 0; i < numStates; i++)            {                bwd[i, t] = 0;                for (int j = 0; j < numStates; j++)                    bwd[i, t] += (bwd[j, t + 1] * TransitionProbabilities[i, j] * EmissionProbabilities[j, o[t + 1]]);            }        }        return bwd;    }