First-order Markov Model

馬可夫模型大意是：選一個狀態作為起點，然後沿著邊隨意走訪任何一個狀態，一直走一直走，沿途累計機率，走累了就停在某個狀態。

熟悉「圖論」的讀者應該很快就能上手，馬可夫模型的外觀看起來就是圖──只不過代數符號多得令人生厭。

建議閱讀：L. Rabiner, A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition, Proceedings of the IEEE, vol. 77, No. 2, February 1989.

State, S

圖片中，每一個圓圈叫做一個「狀態」。馬可夫模型一共有N個狀態，通常標作s₁到s_N。N是我們自行設定的常數。

全部狀態構成的集合，標作大寫S。

Transition Probability, A

每一個狀態都會射出N條邊，這N條邊分別連向每一個狀態，這N條邊的數值都是機率，這N條邊的數值皆介於0到1，這N條邊的數值總和必為1，滿足機率公設。

一條由狀態s_i到狀態s_j的邊，其數值通常標作a_ij。亦可標作條件機率P( s_j | s_i )，意思是現在在狀態s_i、要來去狀態s_j。亦可套用稍後提到的狀態序列，標作P( q_t+1 = j | q_t = i )。

全部邊構成的集合，標作A。通常把A看作一個N×N矩陣。

寫程式時，我們可以利用圖的資料結構adjacency matrix儲存全部的數值。當A有很多數值是零，去掉一些邊之後，也可以改用adjacency lists。

Initial Probability, Π

我們可以選擇任何一個狀態作為起點。機率總和為1。

繪製圖片時，我們可以設計一個虛幻狀態s₀，讓s₀射出N條虛幻邊，分別連向每一個狀態，其數值滿足機率公設。如此一來我們就可以從冥冥之中的s₀開始行動。

這N條虛幻邊的數值通常標作π₁到π_N。亦可標作機率P( s₁ )到P( s_N )。亦可套用稍後提到的狀態序列，標作P( q₁ = 1 )到P( q₁ = N )。

這N條虛幻邊構成的集合，標作Π。通常把Π看作一個N維向量、一個N×1矩陣。

寫程式時，通常我們不會設計一個虛幻的狀態S₀，而是把狀態s₁到s_N重新標作s₀到s_N-1。

const int N = 3; double a[N][N]; double π[N];

State Sequence, q₁ q₂ ...... q_T

選定起點後，可以沿著邊，到處走來走去，一路走過T個狀態以及T-1條邊之後，就停下來。T是我們自行設定的常數。

一條路徑，通常標作q₁ q₂ ...... q_T，為途經的狀態編號。

一條路徑的機率，可以寫成π_q1 × a_q1q2 × ...... × a_qT-1qT，也可以寫成P( s_q1 ) × P( s_q2 | s_q1 ) × ...... × P( s_qT-1 | s_qT )。

const int N = 3; double π[N], a[N][N]; // Markov Model double probability(int* q, int T) { double p = π[q[0]]; for (int i=1; i<T; ++i) p *= a[q[i-1]][q[i]]; return p; }

延伸閱讀：time variant

馬可夫模型的A的數值，會隨時間而改變。

時刻t，一條由狀態s_i到狀態s_j的邊，其數值通常標作a_tij。

平常不會搞到這麼複雜啦。一般只用time invariant。

延伸閱讀：order

馬可夫模型的A的數值，會受到來源狀態影響。

如果只受到上一個狀態影響，稱做first-order，也就是我們慣用的：一條由狀態s_i到狀態s_j的邊，其數值標作a_ij，亦可標作P( s_j | s_i )，亦可標作P( q_t+1 = j | q_t = i )。

如果受到上一個狀態、上上個狀態影響，稱做second-order：其數值標作a_ijk，亦可標作P( s_k | s_i s_j )，亦可標作P( q_t+1 = k | q_t = j , q_t-1 = i )。

如果受到以前所有狀態影響，我不知道要稱做甚麼：可標示成條件機率P( q_t+1 = j | q_t = i_t , q_t-1 = i_t-1 , ...... , q₀ = i₀ )。高中數學提到的Markov Model其實是這個版本。

平常不會搞到這麼複雜啦。一般只用first-order。

Observation, V

隱藏馬可夫模型添加了一個新要素：每當造訪一個狀態，就立刻從M個值當中，噴出其中一個值。每一個狀態都是噴出相同的M種值，這M個值通常標作v₁到v_M。M是我們自行設定的常數。

全部觀察構成的集合，標作大寫V。

Observation Probability, B

每一個狀態噴出這M種值的機率都不相同。

狀態s_i噴出v_k的機率，通常標作b_i( k )或者簡單標作b_ik。亦可標作條件機率P( v_k | s_i )，意思是現在在狀態s_i、噴出觀察v_k。亦可套用狀態序列與觀察序列，標作P( o_t = k | q_t = i )。

每個狀態各自的噴出機率構成的集合，標作B。通常把B看作N個函數，或者簡單地看作一個N×M矩陣。

Observation Sequence, o₁ o₂ ...... o_T

走T步後，一路上T個狀態個別噴出的值的編號。

有了狀態序列與觀察序列，一條路徑的機率，可以寫成π_q1 × b_q1( o₁ ) × a_q1q2 × b_q2( o₂ ) × ...... × a_qT-1qT × b_qT( o_T )，也可以寫成P( s_q1 ) × P( v_o1 | s_q1 ) × P( s_q2 | s_q1 ) × P( v_o2 | s_q2 ) × ...... × P( s_qT-1 | s_qT ) × P( v_oT | s_qT )。我想各位差不多眼花撩亂了。名可名，非常名，若能領會原理就不用刻意背誦代數了。

const int N = 3, M = 3, T = 15; double π[N], a[N][N], b[N][M]; // HMM double probability(int* q, int* o, int T) { double p = π[q[0]] * b[q[0]][o[0]]; for (int i=1; i<T; ++i) p *= a[q[i-1]][q[i]] * b[q[i]][o[i]]; return p; }

隱藏馬可夫模型的特色就是：我們只看到了觀察序列（果），但是我們看不到狀態序列（因）；我們只看到了依序噴出的T個值，但是我們看不到一路走過的是哪T個狀態。

「隱藏」二字便是指行蹤被隱藏了，狀態序列被隱藏了。

接下來要討論隱藏馬可夫模型的三個基本問題，以及演算法。

1. Evaluation Problem: Forward-backward Algorithm

看到一個觀察序列o₁ o₂ ...... o_T，但是看不到狀態序列s₁ s₂ ...... s_T的情況下，找出所有可能的路徑的機率的總和。

對於一個觀察序列來說，狀態序列有各式各樣的可能性，一共有N^T種可能性。

運用窮舉法，時間複雜度為O(N^T * T)。運用「動態規劃」，時間複雜度降低為O(N² * T)。

原理是結合律，比如x×b×c + y×b×c + z×b×c = (x+y+z)×b×c，能加的先加一加。左端結合是forward，右端結合是backward；使用其中一種即可，計算結果都一樣。

forward:
α1  (j) = πj × bj(o1)
αt+1(j) = Σi=1~N[ αt(i) × aij ] × bj(ot+1)

backward:
βT(i) = 1
βt(i) = Σj=1~N[ aij × bj(ot+1) × βt+1(j) ]

const int N = 3, M = 3, T = 15; double π[N], a[N][N], b[N][M]; // HMM double α[T][N]; // 可以簡化成α[2][N] double β[T][N]; // 可以簡化成β[2][N] double forward(int* o, int T) { for (int t=0; t<T; ++t) for (int j=0; j<N; ++j) if (t == 0) α[t][j] = π[j] * b[j][o[t]]; else { double p = 0; for (int i=0; i<N; ++i) p += α[t-1][i] * a[i][j]; α[t][j] = p * b[j][o[t]]; } double p = 0; for (int i=0; i<N; ++i) p += α[T-1][i]; return p; } double backward(int* o, int T) { for (int t=T-1; t>=0; --t) for (int i=0; i<N; ++i) if (t == T-1) β[t][i] = 1.0; else { double p = 0; for (int j=0; j<N; ++j) p += a[i][j] * b[j][o[t+1]] * β[t+1][j]; β[t][i] = p; } double p = 0; for (int j=0; j<N; ++j) p += π[j] * b[j][o[0]] * β[0][j]; return p; }

2. Decoding Problem: Viterbi Algorithm

看到一個觀察序列o₁ o₂ ...... o_T，但是看不到狀態序列s₁ s₂ ...... s_T的情況下，從所有可能的路徑當中，找出機率最大的一條路徑，以及其機率。

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

forward:
δ1  (j) = πj × bj(o1)
δt+1(j) = maxi=1~N[ δt(i) × aij ] × bj(ot+1)

path tracing:
ψt(j) = arg maxi=1~N[ δt-1(i) × aij ] × bj(ot)
      = arg maxi=1~N[ δt-1(i) × aij ]

const int N = 3, M = 3, T = 15; double π[N], a[N][N], b[N][M]; // HMM double δ[T][N]; // 可以簡化成δ[2][N] int ψ[T][N]; double decode(int* o, int T, int* q) { for (int t=0; t<T; ++t) for (int j=0; j<N; ++j) if (t == 0) δ[t][j] = π[j] * b[j][o[t]]; else { double p = -1e9; for (int i=0; i<N; ++i) { double w = δ[t-1][i] * a[i][j]; if (w > v) p = w, ψ[t][j] = i; } δ[t][j] = p * b[j][o[t]]; } double p = -1e9; for (int j=0; j<N; ++j) if (δ[T-1][j] > v) p = δ[T-1][j], q[T-1] = j; for (int t=T-1; t>0; --t) q[t-1] = ψ[t][q[t]]; return p; }

3. Learning Problem: EM Algorithm

給定一個觀察序列o₁ o₂ ...... o_T，更新ABΠ使得Evaluation Problem算得的機率盡量大。

更新的原理，採用了「Maximum Likelihood Estimation」，以樣本平均值作為分布平均值，出現這些樣本的機率就是最大。

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

              αt(i) × βt(i)
γt(i) = ————————————————————————— = Σj=1~N ξt(i,j)
         Σi=1~N [ αt(i) × βt(i) ]

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

                  αt(i) × aij × bj(ot+1) × βt+1(j)
ξt(i,j) = —————————————————————————————————————————————————
           Σi=1~N Σj=1~N [ αt(i) × aij × bj(ot+1) × βt+1(j) ]

更新：

πi    = γ1(i)
aij   = Σt=1~T-1 [ ξt(i,j) ] ÷ Σt=1~T-1 [ γt(i) ]
bj(k) = Σt=1~T, ot=k [ γt(i) ] ÷ Σt=1~T [ γt(i) ]

const int N = 3, M = 3, T = 15; double π[N], a[N][N], b[N][M]; // HMM double α[T][N], β[T][N]; // evalution problem double δ[T][N]; int ψ[T][N]; // decoding problem double γ[T][N], ξ[T][N][N]; // learning problem void learn(int* o, int T) { forward(o, T); backward(o, T); for (int t=0; t<T; ++t) { double p = 0; for (int i=0; i<N; ++i) p += α[t][i] * β[t][i]; assert(p != 0); for (int i=0; i<N; ++i) γ[t][i] = α[t][i] * β[t][i] / p; } for (int t=0; t<T-1; ++t) { double p = 0; for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) p += α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j]; assert(p != 0); for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) ξ[t][i][j] = α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j] / p; } // 更新Π for (int i=0; i<N; ++i) π[i] = γ[0][i]; // 更新A for (int i=0; i<N; ++i) { double p2 = 0; for (int t=0; t<T-1; ++t) p2 += γ[t][i]; assert(p2 != 0); for (int j=0; j<N; ++j) { double p1 = 0; for (int t=0; t<T-1; ++t) p1 += ξ[t][i][j]; a[i][j] = p1 / p2; } } // 更新B for (int i=0; i<N; ++i) { double p[M] = {0}, p2 = 0; for (int t=0; t<T; ++t) { p[o[t]] += γ[t][i]; p2 += γ[t][i]; } assert(p2 != 0); for (int k=0; k<M; ++k) b[i][k] = p[k] / p2; } }

延伸閱讀：取log

寫程式時，機率都是小於一的數字，連乘之後數字越來越小。然而，計算機的浮點數，精確度是有極限的，當T很大時，連乘之後那就變成零了。所以實作時我們會取log，連乘也就變成了連加，避免連乘之後變成零的窘境。

舉例來說，一條路徑的機率，取log之後，可以寫成：log( π_q1 × a_q1q2 × b_q1( o₁ ) × ...... ) = log( π_q1 ) + log( a_q1q2 ) + log( b_q1( o₁ ) ) + ......。預先把ABΠ的每個數值取log即可。

取log之後，處理decoding problem沒有什麼大問題，比較麻煩的是evaluation problem與learning problem，除了乘法運算還有加法運算。實數乘法化作了log加法，那麼實數加法怎麼辦呢？可以使用下列公式：

利用log(p)與log(q)，求出log(p+q)，一般讓底數等於10。b可以是任意數。

if p ≥ q
  logb (p + q)  = log p + logb (1 + blogb q - logb p)
else
  logb (p + q)  = log q + logb (1 + blogb p - logb q)

知名的隱藏馬可夫模型套件HTK是這麼實作的：

#define LZERO (-1.0E10) // log(0) #define LSMALL (-0.5E10) // log values < LSMALL are set to LZERO #define minLogExp -log(-LZERO) // ~=-23 double LogAdd(double x, double y) { double temp, diff, z; if (x < y) { temp = x; x = y; y = temp; } diff = y-x; // notice that diff <= 0 if (diff < minLogExp) // if y' is far smaller than x' return (x < LSMALL) ? LZERO : x; else { z = exp(diff); return x + log(1.0 + z); } }

延伸閱讀：smoothing

Learning Problem針對沒出現的觀察值，更新之後機率是零。往後，Evalution Problem與Decoding Problem針對此觀察值，算得機率均是零。也就是說，一旦變成零，就無法恢復成非零了。

更新為零，可以改成更新為一個很小但是不等於零的數字，記得維持機率總和等於一。

只有離散版本HMM有此問題，稍後介紹的連續版本HMM就無此問題。

Pattern Recognition

   兩組調性不同的資料集
   Data Set 1    Data Set 2
1. ABBCABCAABC   BBBCCBC
2. ABCABC        CCBABB
3. ABCAABC       AACCBBB
4. BBABCAB       BBABBAC
5. BCAABCCAB     CCAABAB
6. CACCABCA      BBBCCBAA
7. CABCABCA      ABBBBABA
8. CABCA         CCCCC
9. CABCA         BBAAA

每一個資料集各自建立一個HMM，實施Training Problem訓練旗下所有資料，每筆資料輪流訓練一次，然後循環L次。如果不採輪流訓練，每筆資料都是連著訓練L次的話，HMM熟記了後面幾個，反而淡忘了前面幾個。

const int N = 3, M = 3, T = 15, L = 20; struct HMM { double π[N], a[N][N], b[N][M]; double forward(int* o, int T); double backward(int* o, int T); double decode(int* o, int T, int* q); void learn(int* o, int T); } hmm[2]; // 每個HMM可以共用的記憶體。 double α[2][N], β[2][N]; // evalution problem double δ[2][N]; int ψ[T][N]; // decoding problem double γ[T][N], ξ[T][N][N]; // learning problem void train(string dataset[2][9]) { for (int i=0; i<2; ++i) for (int l=0; l<L; ++l) for (int j=0; j<9; ++j) { // 把資料轉換成離散數值。最好可以預先轉換好。 int o[T]; for (int k=0; k<dataset[i][j].size(); ++k) o[k] = dataset[i][j][k] - 'A'; hmm[i].learn(o, dataset[i][j].size()); } }

如果每筆資料是不定期收到的、訓練好的HMM隨時要用於辨識，就無法輪流訓練。解決方式是記住之前的Training Problem表格，以比重求得這次Training Problem的平均值。

          1 × γ1(i) + n × γ1(n)(i)
πi(n+1) = ——————————————————————————
                   1 + n

           1 × Σt=1~T-1 [ ξt(i,j) ] + n × Σt=1~T-1 [ ξt(n)(i,j) ]
aij(n+1) = ——————————————————————————————————————————————————————
               1 × Σt=1~T [ γt(i) ] + n × Σt=1~T [ γt(n)(i) ]

             1 × Σt=1~T, ot=k [ γt(i) ] + n × Σt=1~T, ot=k [ γt(n)(i) ]
bj(n+1)(k) = ———————————————————————————————————————————————————————
                  1 × Σt=1~T [ γt(i) ] + n × Σt=1~T [ γt(n)(i) ]

Pattern Recognition

訓練好的HMM馬上可用於辨識。要辨識一筆資料，就窮舉每一個HMM，看看該筆資料最符合哪個HMM。所謂最符合，就是實施Evaluation Problem，取機率最大者，作為辨識結果。

為了加速，通常改用Decoding Problem。取log的情況下，Evaluation Problem要處理對數相加，而Decoding Problem不必，因此Decoding Problem效率略勝一籌。

經實驗觀察，採用Evaluation Problem與採用Decoding Problem的辨識結果幾乎相同。至於你信不信，我反正信了。

int recognize(int* o, int T) { int ans = -1; double p = -1e9; for (int i=0; i<2; ++i) { double pp = hmm[i].decode(o, T, q); if (pp > p) p = pp, ans = i; } return ans; }

Guassian Mixture Model

之前說，每個狀態都會噴出M種固定的值。然而現實生活中，觀察值通常是實數，也通常有誤差，不見得恰好就是這M種值。

一個直覺的解決方案是假設誤差呈常態分布：以常態分布的平均值μ，代表正確的噴出數值；以常態分布的變異數σ²，控制誤差範圍。M個觀察值分別製作M個常態分布，標作N( v , μ₁ , σ₁² )到N( v , μ_M , σ_M² )。

狀態不再噴出k種值，而是有k種噴出管道。觀察值不再是固定的、離散的數值v₁到v_M，而是變動的、連續的數值v。

使用第k個噴出管道的機率，仍標作b_i( k )。在第k個噴出管道，噴出數值v的機率，是一個常態分布函數N( v , μ_k , σ_k² )；v代入一個實際數值，就能計算其噴出機率。

當各個常態分布的平均值很接近，或者變異數很大時，不同的噴出管道可能噴出相同數值。狀態s_i噴出數值v的機率，M個噴出管道都必須累計，為Σ_i=1~M [ b_i( k ) × N( v , μ_k , σ_k² ) ]。

綜觀整個過程，是把M個常態分布以不同比重疊加起來，合成一個分布──此概念稱作高斯混合模型，可以應用在許多地方，連續版本的HMM便是一例。另外，通常會將b_i( k )重新標作c_ik，以呼應加權比重的概念。

寫程式時，我們不會真的疊加起來，而是分別紀錄M個高斯混合模型的平均值和變異數。

每個狀態噴出不同的觀察值

套用高斯混合模型，觀察值從原先的離散數值變成了連續數值。換個角度想，套用高斯混合模型，就能製造任何一種我們想要的連續分布。

每個狀態可以使用不同平均值、不同變異數的常態分布，甚至可以使用不同數量的常態分布。狀態s_i使用的常態分布，標作N( v , μ_ik , σ_ik² )。

觀察值推廣到高維度

觀察值可以從一個值推廣為一個向量。其實也就是每個維度分開處理罷了。

至於Learning Problem，則要額外更新常態分布的平均值、變異數、加權比重了：

　　　路徑第一步，穿過狀態si。
πi  =  γ1(i)

　　　分子是穿過邊aij。分母是穿過狀態si。
       Σt=1~T-1 [ ξt(i,j) ]
aij = —————————————————————
         Σt=1~T [ γt(i) ]

　　　分子是穿過狀態si、使用第k個噴出管道。分母是穿過狀態si。
          Σt=1~T [ γtk(i) ]
cik = ————————————————————————
       Σt=1~T Σk=1~M [ γtk(i) ]

　　　分子是穿過狀態si、使用第k個噴出管道、噴出數值的平均值。
       Σt=1~T [ γtk(i) ] * ot
μik = ————————————————————————
          Σt=1~T [ γtk(i) ]

　　　分子是穿過狀態si、使用第k個噴出管道、噴出數值的共變異矩陣。
       Σt=1~T [ γtk(i) ] * (ot - μik)(ot - μik)T
Σik = ——————————————————————————————————————————
                  Σt=1~T [ γtk(i) ]

至此產生了連續版本的HMM，已經可以應付各種情況了！