如果你不想看那麼多字想直接看程式碼,可以直接點上面的連結,程式碼都有註解,應該不難懂。
不過我都寫了還是希望可以看一下,拜託啦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 | t = 1 | t = 2 | t = 3 | t = 4 | |
---|---|---|---|---|---|
state = 0 | |||||
state = 1 |
所以要幹嘛就很清楚了,把這個表格填滿就好了。
接著我們要知道計算
在計算
公式會調整成這樣:
你可能會問說那個scaling coefficient在哪?其實就是
再來我們一個一個步驟來重新解析:
def Forward(O, a, b, pi, isBW=False):
# init parameters
# scaled alpha (hat alpha)
scaled_alpha = np.zeros((O.shape[0], a.shape[0]))
alpha = np.zeros(a.shape[0]) # alpha
C = np.zeros(O.shape[0]) # scaling factor
C[0] = np.sum(pi * b[:, O[0]])
scaled_alpha[0, :] = (pi * b[:, O[0]]) / C[0] # init scaled alpha
P = 0 # probability of the observation
# compute alpha, t: time t
for t in range(1, O.shape[0]):
for j in range(a.shape[0]):
alpha[j] = scaled_alpha[t - 1] @ a[:, j] * b[j, O[t]]
C[t] += alpha[j]
C[t] = 1 / C[t]
P += math.log2(C[t])
scaled_alpha[t, :] = alpha * C[t] # update alpha
P *= -1
print("Oberservation Probability: ", P)
## if is called by Baum-Welch, return scaled_alpha
## if is called by other function, return P
if isBW:
return scaled_alpha
else:
return P
因為篇幅關係程式碼解說請點這裡
Backward其實想法跟Forward很相似,不過他是從最後一個
為了與Forward的
Backward的終極目標是求得
我們把
t = 0 | t = 1 | t = 2 | t = 3 | t = 4 | |
---|---|---|---|---|---|
state = 0 | |||||
state = 1 |
既然
公式會調整成這樣:
拆成一個一個步驟來看:
def Backward(O, a, b):
# init beta
scaled_beta = np.zeros((O.shape[0], a.shape[0])) # scaled beta
beta = np.zeros(a.shape[0])
C = np.zeros(O.shape[0])
P = 0 # probability of the backward observation
# setting beta(T) = 1
scaled_beta[O.shape[0] - 1] = np.ones((a.shape[0]))
# compute beta (from time t-1 to 1)
for t in range(O.shape[0] - 2, -1, -1):
for j in range(a.shape[0]):
beta[j] = (scaled_beta[t + 1] * b[:, O[t + 1]]) @ a[j, :]
C[t + 1] += beta[j]
C[t + 1] = 1 / C[t + 1]
P += math.log2(C[t + 1])
scaled_beta[t, :] = beta * C[t + 1]
P *= -1
# print("Oberservation Backward Probability: ", P)
return scaled_beta
因為篇幅關係程式碼解說請點這裡
講了這麼多終於要進入本篇的重頭戲-BaumWelch演算法了,BW演算法正是利用前面的Forward和Backward所計算出來的
在HMM的概念中,我們根本不會知道最一開始的A和B是多少,所以才需要BW演算法來求得最有可能的A和B。
BW中有兩個很重要的參數,一個是
先來講解
公式:
接著是
公式:
此外,
和 是有關係的: ,意即把當前時間點同個狀態的 全加起來就是
恭喜,我們現在有了新的A和B了,我們只要不斷重複這步驟就可以找到更好的A和B,至於要跑幾次就交給你決定了。
def BaumWelch(O, a, b, pi, n_iter=100):
N = a.shape[0]
T = len(O) # total time
for i in range(n_iter):
alpha = Forward(O, a, b, pi, True)
# print("alpha:", alpha)
beta = Backward(O, a, b)
# print("beta:", beta)
xi = np.zeros((N, N, T - 1)) # init xi
# compute xi
# scaled alpha and beta wont effect xi
# xi: probability of being in state i at time t and state j at time t + 1
for t in range(T - 1):
denominator = np.dot(
np.dot(alpha[t, :].T, a) * b[:, O[t + 1]].T, beta[t + 1, :])
for i in range(N):
numerator = alpha[t, i] * a[i, :] * \
b[:, O[t + 1]].T * beta[t + 1, :].T
xi[i, :, t] = numerator / denominator
# use xi to compute gamma
# scaled alpha and beta wont effect xi
# gamma : probability of being in state i at time t (fix i and t then sum all j)
gamma = np.sum(xi, axis=1)
# use xi and gamma to update a
# sum xi and gamma over time t
# new a = sum(xi) / sum(gamma)
# because sum(xi, 2) is two demension and sum(gamma, axis=1) is one demension, gamma need to be expand to two demension
# fix row, reshape column
a = np.sum(xi, 2) / np.sum(gamma, axis=1).reshape((-1, 1))
# Add additional T'th element in gamma (because xi is T-1)
# hsack: stack arrays in sequence horizontally (column wise)
# use the origin gamma formula to compute the last gamma
denominator = alpha[T - 1] @ beta[T - 1]
gamma = np.hstack((gamma, np.expand_dims(
(alpha[T - 1] * beta[T - 1]) / denominator, axis=1)))
# use gamma to update b
# K is the number of observation types (four in this case)
K = b.shape[1]
denominator = np.sum(gamma, axis=1)
for l in range(K):
b[:, l] = np.sum(gamma[:, O == l], axis=1)
b = np.divide(b, denominator.reshape((-1, 1)))
return a,
同樣因為篇幅關係我把程式碼解說放在這裡
一般在講解HMM的時候,都會從時間點t=1開始,
但python的陣列是從index=0開始,所以我們會從時間點t=0開始。
這裡用
和
程式碼在實作時不會用到這行,因為如果是計算BW演算法只需要