程式扎記

來源自這裡
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 代碼:

view plaincopy to clipboardprint?
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
        for (int j=0; j
            if (t == 0)  
                α[t][j] = π[j] * b[j][o[t]];  
            else  
            {  
                double p = 0;  
                for (int i=0; i
                    p += α[t-1][i] * a[i][j];  
                α[t][j] = p * b[j][o[t]];  
            }  
  
    double p = 0;  
    for (int i=0; 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
            if (t == T-1)  
                β[t][i] = 1.0;  
            else  
            {  
                double p = 0;  
                for (int j=0; j
                    p += a[i][j] * b[j][o[t+1]] * β[t+1][j];  
                β[t][i] = p;  
            }  
  
    double p = 0;  
    for (int j=0; j
        p += π[j] * b[j][o[0]] * β[0][j];  
    return p;  
}  

Decoding Problem: Viterbi Algorithm
看到一個觀察序列 o1 o2 ...... oT ，但是看不到狀態序列 s1 s2 ...... sT 的情況下，從所有可能的路徑當中，找出機率最大的一條路徑，以及其機率。

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

- C 代碼

view plaincopy to clipboardprint?
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
        for (int j=0; j
            if (t == 0)  
                δ[t][j] = π[j] * b[j][o[t]];  
            else  
            {  
                double p = -1e9;  
                for (int i=0; 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
        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;  
}  

Learning Problem: EM Algorithm
給定一個觀察序列 o1 o2 ...... oT ，更新 ABΠ 使得 Evaluation Problem 算得的機率盡量大. 更新的原理，採用了「 Maximum Likelihood Estimation 」，以樣本平均值作為分布平均值，出現這些樣本的機率就是最大:

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

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

更新:

- C 代碼:

view plaincopy to clipboardprint?
const int N = 3, M = 3, T = 15;  
double π[N], a[N][N], b[N][M];  // HMM  
double α[T][N], β[T][N];        // evaluation 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
    {  
        double p = 0;  
        for (int i=0; i
            p += α[t][i] * β[t][i];  
        assert(p != 0);  
  
        for (int i=0; i
            γ[t][i] = α[t][i] * β[t][i] / p;  
    }  
  
    for (int t=0; t1; ++t)  
    {  
        double p = 0;  
        for (int i=0; i
            for (int j=0; j
                p += α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j];  
        assert(p != 0);  
  
        for (int i=0; i
            for (int j=0; j
                ξ[t][i][j] = α[t][i] * a[i][j] * b[j][o[t+1]] * β[t+1][j] / p;  
    }  
  
    // 更新Π  
    for (int i=0; i
        π[i] = γ[0][i];  
  
    // 更新A  
    for (int i=0; i
    {  
        double p2 = 0;  
        for (int t=0; t1; ++t)  
            p2 += γ[t][i];  
        assert(p2 != 0);  
  
        for (int j=0; j
        {  
            double p1 = 0;  
            for (int t=0; t1; ++t)  
                p1 += ξ[t][i][j];  
            a[i][j] = p1 / p2;  
        }  
    }  
  
    // 更新B  
    for (int i=0; i
    {  
        double p[M] = {0}, p2 = 0;  
        for (int t=0; t
        {  
            p[o[t]] += γ[t][i];  
            p2 += γ[t][i];  
        }  
        assert(p2 != 0);  
  
        for (int k=0; k
            b[i][k] = p[k] / p2;  
    }  
}  

延伸閱讀：取 log
寫程式時，機率都是小於一的數字，連乘之後數字越來越小。然而，計算機的浮點數，精確度是有極限的，當 T 很大時，連乘之後那就變成零了。所以實作時我們會取 log ，連乘也就變成了連加，避免連乘之後變成零的窘境. 取 log 之後，處理 decoding problem 沒有什麼大問題，比較麻煩的是 evaluation problem 與 learning problem ，除了乘法運算還有加法運算。實數乘法化作了 log 加法，那麼實數加法怎麼辦呢？可以使用下列公式：

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

view plaincopy to clipboardprint?
#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);  
    }  
}  

Supplement:
* Viterbi algorithm : Finding the most likely sequence of hidden states

程式扎記

標籤

2014年5月11日星期日

[ Algorithm ] Hidden Markov Model

4 則留言:

[Git 常見問題] error: The following untracked working tree files would be overwritten by merge

檢舉濫用情形

學習筆記

標籤

2014年5月11日 星期日