Markov chain model (馬可夫鍊模型) 與求解演算法 forward/backward algorithm

在電腦科學中，graphical model 是一個用來定義機率模型的好工具
而馬可夫鍊模型是graphical mode中常見的一個機率模型

最近看了 Mark Schmidt 所寫的 Matlab 函式庫和教學文件
覺得作者的講解很清楚
只要有基本的機率知識就能夠了解
有興趣的人可以看一看

首先推薦其中一篇說明 Markov chain model 教學文件
以 computer science 學生畢業後的出路為範例
每個人在他畢業後的60年內
每一年都會處於七種狀態之一
[業界, 唸博士, 玩電玩, 業界(有博士學位), 學界(有博士學位), 玩電玩(有博士學位), 離開電腦科學]
請問一個人在畢業後的第n年處於某一個狀態的機率是多少?


以下是我個人的心得
Markov chain 模型中所要描述的為一組隨機變數的機率
令這些隨機變數為 `X=X_1,X_2,...,X_N`
這些隨機變數的 joint probability 可以用一個機率函數 表示
`p(x)=p(x_1)p(x_2|x_1)p(x_3|x_2)...p(x_N|x_{N-1})`

而這個函式可以用更一般化的形式來表示
`p(x)=\frac{1}{Z}f_1(x_1,x_2)f_2(x_2,x_3)...f_{N-1}(x_{N-1},x_N)` (1)
由於加上了常數 Z , 此處的 f 不再限定是機率函數
以上的形式表現出Markov chain的性質:
由於每個隨機變數只和他相鄰的隨機變數相關
所以機率函式可以被分解為多個函數相乘

對於 `p(x)` 我們通常會感興趣的是某一個隨機變數處於某個狀態的機率
比如說 `x_2`=業界 (狀態為1)
以機率的術語來說就是 marginal probability
由於我們只考慮 `x_2` 的狀態
因此要把其他的變數的所有可能性都加總起來
`p(x_2=1)=\sum_{x_1,x_3,x_4,...)\frac{1}{Z}f_1(x_1,x_2)f_2(x_2,x_3)...f_{N-1}(x_{N-1},x_N)`
乍看之下其餘 N-1 個變數有 `2^{N-1}` 種可能
所以要計算這個總和需要非常大量的計算
但實際上卻是很容易計算的
關鍵就在於 

分配律

先看一個簡單的例子
`a_1b_1+a_1b_2+a_2b_1+a_2b_2=(a_1+a_2)(b_1+b_2)`
左邊:先乘後加 -- 需要四個乘法和三個加法
右邊:先加後乘 -- 需要一個乘法和兩個加法
可以看出先加後乘在計算上比較有效率
是當變數增加的時候更為明顯

再回頭看(1)
可以發現是很相似的情況
先做 N-2 個乘法再加起來
除了一點不同:函數 f_i 牽涉到兩個變數
並無法把 `\sum_{x_1,\cdots \x_N}` 拆開
不過 當任何一個變數的值被限定時
就可以拆解開來
以上面的例子來說
令 `x2=1`
則我們可以把乘積分成兩段
`p(x_2=1)=\frac{1}{Z} [\sum_{x_1}f_1(x_1,x_2=1)][\sum_{x_3,...x_N}f_2(x_2=1,x_3)...f_{N-1}(x_{N-1},x_N)]`

到此已經解開最困難的關鍵
我們只需要分別計算以下兩個數字即可
`\sum_{x_1}f_1(x_1,x_2=1)` (3)
`\sum_{x_3,...x_N}f_2(x_2=1,x_3)...f_{N-1}(x_{N-1},x_N)` (4)

其中用來計算(3)的方法便稱之為 forward algorithm 為一種動態規劃的技巧 (dynamic programming)

forward algorithm
現在來計算 (3)
令 `Z_{i,j}` 為 `\sum_{x_1,...\x_{i-1}}f_1(x_1,x_2)f_2(x_2,x_3)...f_{i-1}(x_{i-1},x_i=j)` (3)

如此可以利用遞回的定義來計算
`Z_{i,j}=\sum_k Z_{i-1,k}f_i(k,j)`
`Z_{0,j}` 為常數 需要先定義好

有了`Z_{i,j}`之後
要計算 `p(x_i)` 便很容易
`p(x_i=j) = \frac{Z_{i,j}}{\sum_k Z_{i,k}}`

實務技巧

`Z_{i,j}`的值可能會非常大或非常小 此時難以用一般的變數來表示
比如說浮點述的精確度可能會不足
解決的方法在於加上兩個輔助的變數 `\kappa_i`, `a_{i,j}`
令 `Z_{i,j} = a_{i,j}\prod_{k=1}^i \kappa_k`
`\kappa_i` 的設定方式為使得 `\sum_k a_{i,k}=1`
如此一來 `a_{i,j}` 的範圍必定在 1 以內

此外計算 `p(x_i)` 時 `\kappa_i` 的值可以省去:
`p(x_i=j) = \frac{Z_{i,j}}{\sum_k Z_{i,k}}`
`p(x_i=j) = \frac{a_{i,j}\prod_{k=1}^i \kappa_k}{\sum_k a_{i,k}\prod_{k=1}^i \kappa_k}`
`p(x_i=j) = \frac{a_{i,j}}{\sum_k a_{i,k}}`

backward algorithm
和 forward algorithm 相似 差別只在於是從右邊算過來
在此就不贅述

經過以上的說明 有興趣的人可以看 Mark Schmidt 所寫的 Matlab 函式庫和教學文件
有完整的程式 可以實際動手玩玩看