隱馬爾可夫模型（五）——隱馬爾可夫模型的解碼問題(維特比算法）

HMM解碼問題

      給定一個觀察序列O=O₁O₂...O_T,和模型μ=（A,B,π），如何快速有效地選擇在一定意義下“最優”的狀態序列Q=q₁q₂...q_T，使該狀態最好地解釋觀察序列。

           

      一種想法是求出每個狀態的概率r_t(i)最大(r_t(i)=P(qt=si,O|μ）)，記q'_t(i)=arg_Qmax(r_t(i))，但是這樣做，忽略了狀態之間的關系，很可能兩個狀態之間的概率為0，即a_{q't(i)q't+1(i)}=0,這樣求得的“最優”狀態序列是不合法的。

      為防止狀態之間轉移概率為0（斷續問題），換一種思路，不是求單個狀態求得最大值，而是求得整個狀態序列最大值，即求

                                   Q'= arg_QmaxP(Q|O,μ）

      此時用維特比算法，先定義下維特比變量δt(i):在時間t，HMM沿着一條路徑到達狀態si，並輸出觀察序列O=O₁O₂...O_t的最大概率:

　　      δ_t(i)=max P(q₁q₂...q_t=s_i,O₁O₂...O_t|μ)

           

                            t                      t+1

      上圖中，對於從t時刻三個到 t+1時刻的狀態1，到底取狀態1,2還是3，不是看單獨狀態1,2還是3的概率，而是看在狀態1,2,3各自的維特比變量值乘以相應的狀態轉換概率，從中選出最大值，假設2時最大，那么記下t+1時刻狀態1之前的路徑是t時刻的狀態2，以此類推。

      δ_t(i)的遞歸關系式: δ_t+1(i)=max_j δ_t(j)*a_ji*b_i(O_t+1),為了記憶路徑，定義路徑變量ψ_t(i)，記錄該路徑上的狀態s_i的前一個狀態。

維特比算法

      step1 初始化：

              δ_t(i) = π_i*b_i(O₁), 1≤i≤N

              ψ_t(i) = 0

      step2 歸納計算：

　　　　  δ_t(i)=max_1≤j≤N δ_t-1(j)*a_ji*b_i(O_t),2≤t≤T;1≤i≤N

             記憶路徑   ψ_t(i) = arg [max_1≤j≤Nδ_t-1(j)*a_ji*b_i(O_t)]

      step3 終結:

            Q_T' = arg max_1≤i≤N [δ_T(i)]

            P'(Q_T') = max_1≤i≤N [δ_T(i)]

　　 step4 路徑回溯:

             q_t'=ψ_t+1(q_t+1') , t=T-1,T-2...1

時間復雜度

      計算某時刻的某個狀態的前向變量需要比較前一時刻的N個狀態，此時時間復雜度為O(N),每個時刻有N個狀態，此時時間復雜度為N*O(N)=O(N²),又有T個時刻，所以時間復雜度為T*O(N²)=O(N²T)。

程序例證

             

       step1 初始化：δ₁(1) = 0.2*0.5=0.1 ，δ₁(2) = 0.4*0.4=0.16， δ₁(3) = 0.4*0.7=0.21

       step2 歸納計算：δ₂(1) =max[0.1*0.5,0.16*0.3,0.21*0.2]*0.6

       ...

      step3 終結：最佳路徑是δ₄(1)δ₄(2)δ₄(3)最大的一個對應的狀態

      step4 回溯：從最后一個狀態往回返

程序代碼

 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int main()
{
        float a[3][3] = {{0.5,0.2,0.3},{0.3,0.5,0.2},{0.2,0.3,0.5}};
        float b[3][2] = {{0.5,0.5},{0.4,0.6},{0.7,0.3}};
        float result[4][3];
        int list[4] = {0,1,0,1};
        int max[4][3];
        float tmp;
        //step1:Initialization
        result[0][0] = 0.2*0.5;
        result[0][1] = 0.4*0.4;
        result[0][2] = 0.4*0.7;
        
        int i,j,k, count = 1, max_node;
        float max_v;
        //step2:歸納運算
        for (i=1; i<=3; i++)
        {
            for(j=0; j<=2; j++)
            {
                tmp = result[i-1][0] * a[0][j] * b[j][list[count]];
                max[i][j] = 0;
                for(k=1; k<=2; k++)
                {
                    if(result[i-1][k] * a[k][j] * b[j][list[count]] > tmp)
                    {
                        tmp = result[i-1][k] * a[k][j]* b[j][list[count]];
                        max[i][j] = k;
                    }
                   result[i][j] = tmp;
                }
                max_v = result[3][0];
                max_node = 0;
                for (k=1; k<=2; k++)
                {
                    if(result[3][k] > max_v)
                    {
                        max_v = result[3][k];
                        max_node = k;
                    }
                }
            }
            count += 1;
        }
        //step3:終結
       for (i=0; i<=3; i++)
        {
            for(j=0; j<=2; j++)
            {
                printf("%d %d     %f\n",i+1,j+1,result[i][j]);

            }
        }
        printf("Pmax=%f\n", max_v);
        printf("step4:%d   \n", max_node+1);
        //step4:回溯
        for(k=3; k>=1; k--)
        {
            printf("step%d:%d  \n",k, max[k][max_node]+1);
            max_node = max[k][max_node];
        }
        return 0;
    }

 

運行結果

        

最終的序列是3 2 2 2