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 函式庫和教學文件
有完整的程式 可以實際動手玩玩看

找出文章中的所有引用

最近在用 Word 寫計畫報告，我在引用參考資料時，是使用 [xx] 這種格式。[xx] 是個交互參照，用來連到文章的引用列表。到目前為止一切順利，不過後來發現引用列表中有不少項目是沒有被引用到的，為了節省空間要，要找出沒有被引用的的項目把他們刪除。
但是自己用眼睛看過30頁的文章實在太累，所以寫了一個 Python script 來幫我找出有用到的參考並加以排序。程式如下:

import re

refdict = {}
fin = open('doc.txt')
for line in fin:
    refs = re.findall('\[\d+\]', line)
    for ref in refs:
        num = int(ref[1:-1])
        refdict[num] = 1

used_ref = refdict.keys()
used_ref.sort()
print used_ref

Python 下用來處理 bibtex 的函式庫 -- pybtex

最近在寫報告的時候，需要和同事整合各自的參考 bibliography)，我本身是採用 JabRef 來做書目管理， JabRef 使用的格式為 bibtex。為了檢查大家的參考書目書目有沒有重複，尋找了一下python 的函式庫，找到以下兩套:


1. python-bibtex
2. pybtex

python-bibtex 是書目管理工具，類似 JabRef，不過需要在 Linux 環境才能運作，所以就放棄了。
pybtex 是用來取代 latex 中的 bibtex 軟體，基本上功能和 bibtex 相同。好消息是 pybtex 有提供方便的函式可以直接剖析 bib 檔，用法大致如下:

from pybtex.database.input import bibtex
parser = bibtex.Parser(encoding="UTF-8")
bib_data = parser.parse_file('foo.bib')
print bib_data.entries.keys()
for k in bib_data.entries.keys():
    print bib_data.entries[k].fields['title']

有了 bibtex 這個工具之後，剩下來就只是比較 title 等簡單的工作了。

語音合成 (一)

上個月的時候，朋友E先生發想做中文的語音合成。
經過搜尋，我們發現目前以工研院的引擎的效果最好 (這邊有線上demo)

不過還是感覺不夠完美，同時也沒有可以調整的參數，
我們的理想是能做出語調和感情的調整。

動手之前我們先做了一番研究，一份較易懂的參考資料如[1]
語音合成牽涉的範圍很廣，大致可以分成兩個部分:
1. 自然語言處理 (Natural Language Processing): 讀取文字，分析出每個字的唸法
2. 數位訊號處理 (Digital Signal Processing): 把獨立的字的發音組合成一個句子

我們決定從第二部分下手，
因為第一部分可以先由人來給訂，而且也無關乎發音是否悅耳

而數位訊號處理技術很多，大致上可以分為兩派
第一派被稱之為 synthesis-by-rule，主要由語音學家發展出來。
第二派被稱之為 synthesis-by-concatenation，和前者的差異為，只採用少量的語言學知識，將語音視為聲音的片段，利用技巧來把片段整合起來。
由於不是語音學的專家，我們決定從第二種方法開始。

首先要有語音資料，目前找到能使用的是 gcin 的語音檔 (ogg格式)
http://ftp.twaren.net/local-distfiles/gcin/ogg.tgz


這個檔案包含約一千個字的發音，對應到國語的所有單字 (包含五種聲調)

原本每一個發音檔是放在各自的目錄下，同時目錄的名稱是中文，所以並不好處理，我先用Python script把檔案轉換成漢語拼音，放在同一個目錄下面。

接著為了由於ogg格式比較少有函式庫可以處理，為了方便起見，我先用 Wav2Mp3 這個軟體把 ogg 轉成 wav 格式。

[1] Thierry DUTOIT, "High-quality text-to-speech synthesis: an overview"

(未完待續)

Blocks -- 電腦視覺研究上的MATLAB framework

前陣子在讀ICCV的paper時，有兩位作者引起我很大的興趣，原因不在於論文的好壞，而是他們把自己做實驗時的研究方法公開在網頁上。

為什麼說這是一件重要的事情呢? 這要先岔題一下:
在電腦視覺的領域，一篇好的研究，除了理論完整、實驗漂亮之外，能否讓其他人重複得到該作者的實驗結果也是非常重要的。

目前在這個領域，已經有公用的測試資料庫供大家使用。想像一下，說假設我提出一套方法，可以找出照片中的行人，要評估我的方法好壞，可以用一組前人已經做過實驗的照片，用他們的結果來比較。而這樣的比較結果，自然是放在我的論文之中。不過問題來了，當有人質疑我的實驗結果時，他們必須用想辦法重現我的實驗結果，不過這是一件很困難的事情。
困難的理由在於，首先你要把我的方法寫成程式，這是一件很花功夫的事，然而當你寫好程式之後還有一個難關，在電腦視覺這個領域中 (電腦科學的很多領域也是)，
一個演算法之中會有很多參數，而參數的設定不同就會造成實驗結果的差異
類比一下，好比說我給你一份食譜，裡面說要加牛肉幾兩、油、鹽、糖幾匙等等
如果我今天漏掉了其中一樣東西的分量，那你就只能自己嘗試，找出最好的分量，
不然味道就會差一點
令人驚訝的是，論文中就是會把很多參數給漏掉
這其實也不能責怪作者，當你的系統大到一個程度時，你自然不可能巨細靡遺的把每個參數寫在論文上面，尤其是當這些參數是跟你的演算法無關的部分。

面對這種情況，
許多論文作者並不會主動提供這些資訊，比如說放在個人網頁上面等等
有時候當你寫信去問時，得到的回覆可能是"抱歉，我的程式不知道放在哪個硬碟了"
這通常不是作者敷衍你，而是一個常見的情形，因為研究題目可能兩三年就會有變動，之前寫的程式往往不能再用(啊哈!你想到軟體的重複利用，但是研究人員不一定具有好的軟體工程能力)

比較好心的作者會把他的程式碼公開，讓其他的研究者可以輕鬆的重複實驗，
這樣子通常可以獲得回報:會提高其他研究者引用他的論文的意願，
而論文的被引用次數，是一個研究者最直接也是最重要的指標(甚至比發表的論文數目還重要)

回到一開始，Brian Fulkerson 和 Andrea Vedaldi 這兩位研究員，他們不只公開程式碼，同時也整理了整個實驗的架構，告訴你可以修改那些地方，讓你能在他們的基礎上來改進這個系統。這才真正的站在巨人的肩膀上! (當然這是指在學術界啦，在業界的話這可是公司賴以生存的寶物，怎能輕易公開!)

(未完待續)