程式扎記: [ Algorithm ] Hidden Markov Model

標籤

2014年5月11日 星期日

[ Algorithm ] Hidden Markov Model

來源自 這裡
Preface:
馬可夫模型大意是:選一個狀態作為起點,然後沿著邊隨意走訪任何一個狀態,一直走一直走,沿途累計機率,走累了就停在某個狀態。熟悉「圖論」的讀者應該很快就能上手,馬可夫模型的外觀看起來就是圖──只不過代數符號多得令人生厭:

建議閱讀: L. Rabiner, A Tutorial on Hidden Markov Models and...vol. 77, No. 2, February 1989.

隱藏馬可夫模型添加了一個新要素:每當造訪一個狀態,就立刻從 M 個值當中,噴出其中一個值。每一個狀態都是噴出相同的 M 種值,這 M 個值通常標作 v1 到 vM 。 M 是我們自行設定的常數。全部觀察構成的集合,標作大寫 V :


Observation Probability, B:
每一個狀態噴出這 M 種值的機率都不相同:


狀態 si 噴出 vk 的機率,通常標作 bi( k ) 或者簡單標作 bik 。亦可標作條件機率 P( vk | si ) ,意思是現在在狀態 si 、噴出觀察 vk 。亦可套用狀態序列與觀察序列,標作 P( ot = k | qt = i ) 。每個狀態各自的噴出機率構成的集合,標作 B 。通常把 B 看作 N 個函數,或者簡單地看作一個 N×M 矩陣.

Hidden Markov Model:
隱藏馬可夫模型 (HMM) 的特色就是:我們只看到了觀察序列(果),但是我們看不到狀態序列(因);我們只看到了依序噴出的 T 個值,但是我們看不到一路走過的是哪 T 個狀態:


「隱藏」二字便是指行蹤被隱藏了,狀態序列被隱藏了。接下來要討論隱藏馬可夫模型的三個基本問題,以及演算法。

Evaluation Problem: Forward-backward Algorithm
看到一個觀察序列 o1 o2 ...... oT ,但是看不到狀態序列 s1 s2 ...... sT 的情況下,找出所有可能的路徑的機率的總和:


對於一個觀察序列來說,狀態序列有各式各樣的可能性,一共有 N^T 種可能性. 運用窮舉法,時間複雜度為 O(N^T * T) 。運用「動態規劃」,時間複雜度降低為 O(N^2 * T) . 原理是結合律,比如x×b×c + y×b×c + z×b×c = (x+y+z)×b×c ,能加的先加一加。左端結合是 forward ,右端結合是 backward ;使用其中一種即可,計算結果都一樣:


- C 代碼:
  1. const int N = 3, M = 3, T = 15;  
  2. double π[N], a[N][N], b[N][M];  // HMM  
  3. double α[T][N]; // 可以簡化成α[2][N]  
  4. double β[T][N]; // 可以簡化成β[2][N]  
  5.   
  6. double forward(int* o, int T)  
  7. {  
  8.     for (int t=0; t
  9.         for (int j=0; j
  10.             if (t == 0)  
  11.                 α[t][j] = π[j] * b[j][o[t]];  
  12.             else  
  13.             {  
  14.                 double p = 0;  
  15.                 for (int i=0; i
  16.                     p += α[t-1][i] * a[i][j];  
  17.                 α[t][j] = p * b[j][o[t]];  
  18.             }  
  19.   
  20.     double p = 0;  
  21.     for (int i=0; i
  22.         p += α[T-1][i];  
  23.     return p;  
  24. }  
  25.   
  26. double backward(int* o, int T)  
  27. {  
  28.     for (int t=T-1; t>=0; --t)  
  29.         for (int i=0; i
  30.             if (t == T-1)  
  31.                 β[t][i] = 1.0;  
  32.             else  
  33.             {  
  34.                 double p = 0;  
  35.                 for (int j=0; j
  36.                     p += a[i][j] * b[j][o[t+1]] * β[t+1][j];  
  37.                 β[t][i] = p;  
  38.             }  
  39.   
  40.     double p = 0;  
  41.     for (int j=0; j
  42.         p += π[j] * b[j][o[0]] * β[0][j];  
  43.     return p;  
  44. }  
Decoding Problem: Viterbi Algorithm
看到一個觀察序列 o1 o2 ...... oT ,但是看不到狀態序列 s1 s2 ...... sT 的情況下,從所有可能的路徑當中,找出機率最大的一條路徑,以及其機率。


這跟上一個問題如出一轍,運用「動態規劃」就可以解決。唯一的差別就是把 Σ 換成 max 而已:


- C 代碼
  1. const int N = 3, M = 3, T = 15;  
  2. double π[N], a[N][N], b[N][M];  // HMM  
  3. double δ[T][N]; // 可以簡化成δ[2][N]  
  4. int ψ[T][N];  
  5.   
  6. double decode(int* o, int T, int* q)  
  7. {  
  8.     for (int t=0; t
  9.         for (int j=0; j
  10.             if (t == 0)  
  11.                 δ[t][j] = π[j] * b[j][o[t]];  
  12.             else  
  13.             {  
  14.                 double p = -1e9;  
  15.                 for (int i=0; i
  16.                 {  
  17.                     double w = δ[t-1][i] * a[i][j];  
  18.                     if (w > v) p = w, ψ[t][j] = i;  
  19.                 }  
  20.                 δ[t][j] = p * b[j][o[t]];  
  21.             }  
  22.   
  23.     double p = -1e9;  
  24.     for (int j=0; j
  25.         if (δ[T-1][j] > v)  
  26.             p = δ[T-1][j], q[T-1] = j;  
  27.   
  28.     for (int t=T-1; t>0; --t)  
  29.         q[t-1] = ψ[t][q[t]];  
  30.   
  31.     return p;  
  32. }  
Learning Problem: EM Algorithm
給定一個觀察序列 o1 o2 ...... oT ,更新 ABΠ 使得 Evaluation Problem 算得的機率盡量大. 更新的原理,採用了「 Maximum Likelihood Estimation 」,以樣本平均值作為分布平均值,出現這些樣本的機率就是最大:


分子是穿越狀態 si 的所有路徑的機率的總和:


分子是穿越邊 aij 的所有路徑的機率的總和:


更新:


- C 代碼:
  1. const int N = 3, M = 3, T = 15;  
  2. double π[N], a[N][N], b[N][M];  // HMM  
  3. double α[T][N], β[T][N];        // evaluation problem  
  4. double δ[T][N]; int ψ[T][N];    // decoding problem  
  5. double γ[T][N], ξ[T][N][N];     // learning problem  
  6.   
  7. void learn(int* o, int T)  
  8. {  
  9.     forward(o, T);  
  10.     backward(o, T);  
  11.   
  12.     for (int t=0; t
  13.     {  
  14.         double p = 0;  
  15.         for (int i=0; i
  16.             p += α[t][i] * β[t][i];  
  17.         assert(p != 0);  
  18.   
  19.         for (int i=0; i
  20.             γ[t][i] = α[t][i] * β[t][i] / p;  
  21.     }  
  22.   
  23.     for (int t=0; t1; ++t)  
  24.     {  
  25.         double p = 0;  
  26.         for (int i=0; i
  27.             for (int j=0; j
  28.                 p += α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j];  
  29.         assert(p != 0);  
  30.   
  31.         for (int i=0; i
  32.             for (int j=0; j
  33.                 ξ[t][i][j] = α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j] / p;  
  34.     }  
  35.   
  36.     // 更新Π  
  37.     for (int i=0; i
  38.         π[i] = γ[0][i];  
  39.   
  40.     // 更新A  
  41.     for (int i=0; i
  42.     {  
  43.         double p2 = 0;  
  44.         for (int t=0; t1; ++t)  
  45.             p2 += γ[t][i];  
  46.         assert(p2 != 0);  
  47.   
  48.         for (int j=0; j
  49.         {  
  50.             double p1 = 0;  
  51.             for (int t=0; t1; ++t)  
  52.                 p1 += ξ[t][i][j];  
  53.             a[i][j] = p1 / p2;  
  54.         }  
  55.     }  
  56.   
  57.     // 更新B  
  58.     for (int i=0; i
  59.     {  
  60.         double p[M] = {0}, p2 = 0;  
  61.         for (int t=0; t
  62.         {  
  63.             p[o[t]] += γ[t][i];  
  64.             p2 += γ[t][i];  
  65.         }  
  66.         assert(p2 != 0);  
  67.   
  68.         for (int k=0; k
  69.             b[i][k] = p[k] / p2;  
  70.     }  
  71. }  
延伸閱讀:取 log
寫程式時,機率都是小於一的數字,連乘之後數字越來越小。然而,計算機的浮點數,精確度是有極限的,當 T 很大時,連乘之後那就變成零了。所以實作時我們會取 log ,連乘也就變成了連加,避免連乘之後變成零的窘境. 取 log 之後,處理 decoding problem 沒有什麼大問題,比較麻煩的是 evaluation problem 與 learning problem ,除了乘法運算還有加法運算。實數乘法化作了 log 加法,那麼實數加法怎麼辦呢?可以使用下列公式:


知名的隱藏馬可夫模型套件 HTK 是這麼實作的:
  1. #define LZERO  (-1.0E10) // log(0)  
  2. #define LSMALL (-0.5E10) // log values < LSMALL are set to LZERO  
  3. #define minLogExp -log(-LZERO) // ~=-23  
  4.   
  5. double LogAdd(double x, double y)  
  6. {  
  7.     double temp, diff, z;  
  8.     if (x < y)  
  9.     {  
  10.         temp = x; x = y; y = temp;  
  11.     }  
  12.     diff = y-x; // notice that diff <= 0  
  13.     if (diff < minLogExp)   // if y' is far smaller than x'  
  14.         return (x < LSMALL) ? LZERO : x;  
  15.     else  
  16.     {  
  17.         z = exp(diff);  
  18.         return x + log(1.0 + z);  
  19.     }  
  20. }  
Supplement:
Viterbi algorithm : Finding the most likely sequence of hidden states


沒有留言:

張貼留言

網誌存檔