如果你不想看那麼多字想直接看程式碼,可以直接點上面的連結,程式碼都有註解,應該不難懂。
不過我都寫了還是希望可以看一下,拜託啦XD
雖說標題是HMM,但這篇會著重在HMM中的Baum-Welch演算法,也就是用來訓練HMM的演算法。這篇文章會用python來實作,並且會用到numpy,所以如果你不熟悉python或numpy,可以先看看這篇文章:Python Numpy Tutorial。
所謂的HMM指的是在不同狀態(status)下,會產生的結果(或是序列),所以很常被使用在基因體序列的預測。例如:在基因體中,有A、T、C、G四種碱基,而在不同的狀態下,會產生不同的碱基,所以我們可以用HMM來預測基因體的碱基序列。
關於HMM的基本概念,可以參考一系列影片,我自己覺得說明的很清楚:Hidden Markov Model。
這裡的重點是會使用到以下幾個HMM的參數:
這裡的轉換機率和發射機率都是一個二維矩陣,依照狀態和序列種類來決定,而初始機率和觀測到的序列都是一個一維陣列。
在介紹Baum-Welch演算法之前我必須先介紹Forward Algorithm,因為Baum-Welch演算法是基於Forward Algorithm來做的。Forward Algorithm是用來計算在某個時間點t,觀測到某個序列O的機率,也就是P(O|λ)。這裡的λ是HMM的參數,包含初始機率、轉換機率、發射機率。
Forward Algorithm的終極目標是求得,也就是P(O|λ)。
要求,首先要先求得,也就是在時間點t=0[註1],觀測到的序列為i的機率,也就是。寫成這樣好像有點難懂,我們畫個表格來看看(這裡假設觀測到的序列O長度為5,所以時間長度也是5,state只有兩種):
t = 0 | t = 1 | t = 2 | t = 3 | t = 4 | |
---|---|---|---|---|---|
state = 0 | |||||
state = 1 |
所以要幹嘛就很清楚了,把這個表格填滿就好了。
接著我們要知道計算的公式
在計算的時候,會有一個問題,就是數字會太小,所以會有underflow的問題(電腦存不到那麼小的數字)。為了解決這個問題,我們可以把乘上一個常數,這樣就不會underflow了,而這個常數稱為scaling coefficient
公式會調整成這樣:
你可能會問說那個scaling coefficient在哪?其實就是這個值,也就是的分母,我們用C來表示。
再來我們一個一個步驟來重新解析:
因為篇幅關係程式碼解說請點這裡
Backward其實想法跟Forward很相似,不過他是從最後一個算回來,這也是他們分別被稱為Forward與Backward的原因。
為了與Forward的做區別,我們把Backward要求的稱為。
Backward的終極目標是求得
我們把搬過來照用。
t = 0 | t = 1 | t = 2 | t = 3 | t = 4 | |
---|---|---|---|---|---|
state = 0 | |||||
state = 1 |
:
既然會遇到underflow的問題,自然也會,所以我們也要使用相同技巧來計算
公式會調整成這樣:
的分母同樣是我們的scaling factor C。
拆成一個一個步驟來看:
因為篇幅關係程式碼解說請點這裡
講了這麼多終於要進入本篇的重頭戲-BaumWelch演算法了,BW演算法正是利用前面的Forward和Backward所計算出來的和來計算更新我們的transition probability和emission probability。
在HMM的概念中,我們根本不會知道最一開始的A和B是多少,所以才需要BW演算法來求得最有可能的A和B。
BW中有兩個很重要的參數,一個是(發音為xi),另一個是,我們正是利用這兩個參數來更新A和B
先來講解,它代表了當time = t時,在state i,⽽當time = t + 1時,在state j的機率(會考慮下個時間點的狀態,所以當前時間點是狀態0的話會產生兩個,分別是狀態0和狀態1,所以在儲存時會用到三維矩陣)
公式:
接著是,也就是當time = t時,在state i的機率(和不同不考慮下一個時間點的機率)。
公式:
此外,和是有關係的: ,意即把當前時間點同個狀態的全加起來就是
恭喜,我們現在有了新的A和B了,我們只要不斷重複這步驟就可以找到更好的A和B,至於要跑幾次就交給你決定了。
同樣因為篇幅關係我把程式碼解說放在這裡
一般在講解HMM的時候,都會從時間點t=1開始,
但python的陣列是從index=0開始,所以我們會從時間點t=0開始。
這裡用來表示機率也是因為underflow的關係。
和不同,沒有初始機率,所以一律將最後的狀態設為1。
程式碼在實作時不會用到這行,因為如果是計算BW演算法只需要矩陣,而單純計算序列機率只要Forward即可。