python 隐马尔科夫_机器学习算法之——隐马尔可夫(Hidden Markov ModelsHMM)原理及Python实现...
前言
上星期寫了Kaggle競(jìng)賽的詳細(xì)介紹及入門指導(dǎo)，但對(duì)于真正想要玩這個(gè)競(jìng)賽的伙伴，機(jī)器學(xué)習(xí)中的相關(guān)算法是必不可少的，即使是你不想獲得名次和獎(jiǎng)牌。那么，從本周開(kāi)始，我將介紹在Kaggle比賽中的最基本的也是運(yùn)用最廣的機(jī)器學(xué)習(xí)算法，很多項(xiàng)目用這些基本的模型就能解決基礎(chǔ)問(wèn)題了。
1 概述
隱馬爾可夫模型(Hidden Markov Model，HMM)是結(jié)構(gòu)最簡(jiǎn)單的動(dòng)態(tài)貝葉斯網(wǎng)，這是一種著名的有向圖模型，主要用于時(shí)序數(shù)據(jù)建模(語(yǔ)音識(shí)別、自然語(yǔ)言處理等)。
假設(shè)有三個(gè)不同的骰子(6面、4面、8面)，每次先從三個(gè)骰子里選一個(gè)，每個(gè)骰子選中的概率為1/3，如下圖所示，重復(fù)上述過(guò)程，得到一串?dāng)?shù)字[1 6 3 5 2 7]。這些可觀測(cè)變量組成可觀測(cè)狀態(tài)鏈。
同時(shí)，在隱馬爾可夫模型中還有一條由隱變量組成的隱含狀態(tài)鏈，在本例中即骰子的序列。比如得到這串?dāng)?shù)字骰子的序列可能為[D6 D8 D8 D6 D4 D8]。
隱馬爾可夫模型示意圖如下所示：
圖中，箭頭表示變量之間的依賴關(guān)系。在任意時(shí)刻，觀測(cè)變量(骰子點(diǎn)數(shù))僅依賴于狀態(tài)變量(哪類骰子)，“觀測(cè)獨(dú)立性假設(shè)”。
同時(shí)，t時(shí)刻的狀態(tài)qt僅依賴于t-1時(shí)刻的狀態(tài)qt-1。這就是馬爾可夫鏈，即系統(tǒng)的下一時(shí)刻的狀態(tài)僅由當(dāng)前狀態(tài)決定不依賴以往的任何狀態(tài)(無(wú)記憶性)，“齊次馬爾可夫性假設(shè)”。
2 隱馬爾可夫模型三要素
對(duì)于一個(gè)隱馬爾可夫模型，它的所有N個(gè)可能的狀態(tài)的集合Q={q1,q2,…,qN}，所有M個(gè)可能的觀測(cè)集合V={v1,v2,…,vM}
隱馬爾可夫模型三要素：狀態(tài)轉(zhuǎn)移概率矩陣A，A=[aij]NxN，aij表示任意時(shí)刻t狀態(tài)qi條件下，下一時(shí)刻t+1狀態(tài)為qj的概率;
觀測(cè)概率矩陣B，B=[bij]NxM，bij表示任意時(shí)刻t狀態(tài)qi條件下，生成觀測(cè)值vj的概率;
初始狀態(tài)概率向量Π，Π=(π1,π2,…,πN),πi表示初始時(shí)刻t=1，狀態(tài)為qi的概率。
一個(gè)隱馬爾可夫模型可由λ=(A, B, Π)來(lái)指代。
3 隱馬爾可夫模型的三個(gè)基本問(wèn)題(1) 給定模型λ=(A, B, Π)，計(jì)算其產(chǎn)生觀測(cè)序列O={o1,o2,…,oT}的概率P(O|λ)；
計(jì)算擲出點(diǎn)數(shù)163527的概率
(2) 給定模型λ=(A, B, Π)和觀測(cè)序列O={o1,o2,…,oT}，推斷能夠最大概率產(chǎn)生此觀測(cè)序列的狀態(tài)序列I={i1,i2,…,iT}，即使P(I|O)最大的I；
推斷擲出點(diǎn)數(shù)163527的骰子種類
(3) 給定觀測(cè)序列O={o1,o2,…,oT}，估計(jì)模型λ=(A, B, Π)的參數(shù)，使P(O|λ)最大；
已知骰子有幾種，不知道骰子的種類，根據(jù)多次擲出骰子的結(jié)果，反推出骰子的種類
這三個(gè)基本問(wèn)題在現(xiàn)實(shí)應(yīng)用中非常重要，例如根據(jù)觀測(cè)序列O={o1,o2,…,oT-1,}推測(cè)當(dāng)前時(shí)刻最有可能出現(xiàn)的觀測(cè)值OT，這就轉(zhuǎn)換成基本問(wèn)題(1)；
在語(yǔ)音識(shí)別中，觀測(cè)值為語(yǔ)音信號(hào)，隱藏狀態(tài)為文字，根據(jù)觀測(cè)信號(hào)推斷最有可能的狀態(tài)序列，即基本問(wèn)題(2)；
在大多數(shù)應(yīng)用中，人工指定參數(shù)模型已變得越來(lái)越不可行，如何根據(jù)訓(xùn)練樣本學(xué)得最優(yōu)參數(shù)模型，就是基本問(wèn)題(3)。
4 三個(gè)基本問(wèn)題的解法
基于兩個(gè)條件獨(dú)立假設(shè)，隱馬爾可夫模型的這三個(gè)基本問(wèn)題均能被高效求解。
4.1 基本問(wèn)題(1)解法
4.1.1 直接計(jì)算法(概念上可行，計(jì)算上不可行)
通過(guò)列舉所有可能的長(zhǎng)度為T的狀態(tài)序列I=(i1,i2,…,iT)，求各個(gè)狀態(tài)序列I與觀測(cè)序列O同時(shí)出現(xiàn)的聯(lián)合概率P(I,O|λ)，然后對(duì)所有可能求和得到P(O|λ)。
狀態(tài)序列I=(i1,i2,…,iT)的概率是P(I|λ)=π1a12a23…a(T-1)T
對(duì)于固定狀態(tài)序列 I，觀測(cè)序O=(o1,o2,…,oT)的概率P(O|I,λ)= b11b22…bTT
I 和 O同時(shí)出現(xiàn)的聯(lián)合概率P(I,O|λ)= P(I|λ) P(O|I,λ)
然后對(duì)所有可能的 I 求和，得到P(O|λ)
這種方法計(jì)算量很大，算法不可行。
4.1.2 前向算法(forward algorithm)(t=1,一步一步向前計(jì)算)
前向概率αt(i)= P(o1,o2,…,ot,ii=qi|λ)，表示模型λ，時(shí)刻 t，觀測(cè)序列為o1,o2,…,ot且狀態(tài)為qi的概率。
(1) 初始化前向概率
狀態(tài)為qi和觀測(cè)值為o1的聯(lián)合概率 α1(i)=π1bi1
(2) 遞推t=1,2,…,T-1
根據(jù)下圖，得到αt+1(i)= [Σj=1,…,Nαt(j)αji]bi(t+1)
(3) 終止 P(O|λ) = Σi=1,…,NαT(i)
前向算法高效的關(guān)鍵是其局部計(jì)算前向概率，根據(jù)路徑結(jié)構(gòu)，如下圖所示，每次計(jì)算直接利用前一時(shí)刻計(jì)算結(jié)果，避免重復(fù)計(jì)算，減少計(jì)算量。
4.1.3 后向算法(backward algorithm)
后向概率βt(i) = P(o1,o2,…,ot,ii=qi|λ)，表示模型λ，時(shí)刻t，從t+1到時(shí)刻T觀測(cè)序列o1,o2,…,ot且狀態(tài)為qi的概率。
(1)初始化后向概率
對(duì)最終時(shí)刻的所有狀態(tài)為qi規(guī)定βT(i) = 1。
(2)遞推t=T-1,T-2,…,1
根據(jù)下圖，計(jì)算t時(shí)刻狀態(tài)為qi，且時(shí)刻t+1之后的觀測(cè)序列為ot+1,ot+2,…，ot的后向概率。只需考慮t+1時(shí)刻所有可能的N個(gè)狀態(tài)qi的轉(zhuǎn)移概率(transition probability) αij,以及在此狀態(tài)下的觀測(cè)概率bi(t+1),然后考慮qj之后的后向概率βt+1(j)，得到
βt(i) = Σj=1,…,Nβt+1(j)αijbi(t+1).
?(3) 終止 P(O|λ) = Σi=1,…,Nβ1(i)πibi1
前向算法高效的關(guān)鍵是其局部計(jì)算前向概率，根據(jù)路徑結(jié)構(gòu)，如下圖所示，每次計(jì)算直接利用前一時(shí)刻計(jì)算結(jié)果，避免重復(fù)計(jì)算，減少計(jì)算量。
4.2 基本問(wèn)題(2)解法
4.2.1 近似算法
選擇每一時(shí)刻最有可能出現(xiàn)的狀態(tài)，從而得到一個(gè)狀態(tài)序列。這個(gè)方法計(jì)算簡(jiǎn)單，但是不能保證整個(gè)狀態(tài)序列的出現(xiàn)概率最大。因?yàn)榭赡艹霈F(xiàn)轉(zhuǎn)移概率為0的相鄰狀態(tài)。
4.2.2 Viterbi算法
使用動(dòng)態(tài)規(guī)劃求解概率最大(最優(yōu))路徑。t=1時(shí)刻開(kāi)始，遞推地計(jì)算在時(shí)刻t狀態(tài)為i的各條部分路徑的最大概率，直到計(jì)算到時(shí)刻T，狀態(tài)為i的各條路徑的最大概率，時(shí)刻T的最大概率即為最優(yōu)路徑的概率，最優(yōu)路徑的節(jié)點(diǎn)也同時(shí)得到。
如果還不明白，看一下李航《統(tǒng)計(jì)學(xué)習(xí)方法》的186-187頁(yè)的例題就能明白算法的原理。
狀態(tài)[3 3 3]極為概率最大路徑。
4.3 基本問(wèn)題(3)解法
4.3.1 監(jiān)督學(xué)習(xí)方法
給定S個(gè)長(zhǎng)度相同的(觀測(cè)序列，狀態(tài)序列)作為訓(xùn)練集(O1,I1)，O2,I3)，…, (Os,Is)使用極大似然估計(jì)法來(lái)估計(jì)模型參數(shù)。
轉(zhuǎn)移概率αij 的估計(jì):樣本中t時(shí)刻處于狀態(tài)i，t+1時(shí)刻轉(zhuǎn)移到狀態(tài)j的頻數(shù)為αij，則
觀測(cè)概率bij和初始狀態(tài)概率πi的估計(jì)類似。
4.3.2 Baum-Welch算法
使用EM算法得到模型參數(shù)估計(jì)式
EM算法是常用的估計(jì)參數(shù)隱變量的利器，它是一種迭代方法，基本思想是：
(1) 選擇模型參數(shù)初始值；
(2) (E步)根據(jù)給定的觀測(cè)數(shù)據(jù)和模型參數(shù)，求隱變量的期望；
(3) (M步)根據(jù)已得隱變量期望和觀測(cè)數(shù)據(jù)，對(duì)模型參數(shù)做極大似然估計(jì)，得到新的模型參數(shù)，重復(fù)第二步。
5 Python代碼實(shí)現(xiàn)
5.1 Baum-Welch算法的程序?qū)崿F(xiàn)
這里提供Baum-Welch算法的代碼實(shí)現(xiàn)。
本來(lái)打算試一下用自己寫的HMM跑一下中文分詞，但很可惜，代碼運(yùn)行的比較慢。所以改成 模擬 三角波 以及 正弦波。
# encoding=utf8
?
import numpy as np
import csv
?
class HMM(object):
def __init__(self,N,M):
self.A = np.zeros((N,N)) # 狀態(tài)轉(zhuǎn)移概率矩陣
self.B = np.zeros((N,M)) # 觀測(cè)概率矩陣
self.Pi = np.array([1.0/N]*N) # 初始狀態(tài)概率矩陣
?
self.N = N # 可能的狀態(tài)數(shù)
self.M = M # 可能的觀測(cè)數(shù)
?
def cal_probality(self, O):
self.T = len(O)
self.O = O
?
self.forward()
return sum(self.alpha[self.T-1])
?
def forward(self):
"""前向算法"""
self.alpha = np.zeros((self.T,self.N))
?
# 公式 10.15
for i in range(self.N):
self.alpha[0][i] = self.Pi[i]*self.B[i][self.O[0]]
?
# 公式10.16
for t in range(1,self.T):
for i in range(self.N):
sum = 0
for j in range(self.N):
sum += self.alpha[t-1][j]*self.A[j][i]
self.alpha[t][i] = sum * self.B[i][self.O[t]]
?
def backward(self):
"""后向算法"""
self.beta = np.zeros((self.T,self.N))
?
# 公式10.19
for i in range(self.N):
self.beta[self.T-1][i] = 1
?
# 公式10.20
for t in range(self.T-2,-1,-1):
for i in range(self.N):
for j in range(self.N):
self.beta[t][i] += self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]
?
def cal_gamma(self, i, t):
"""公式 10.24"""
numerator = self.alpha[t][i]*self.beta[t][i]
denominator = 0
?
for j in range(self.N):
denominator += self.alpha[t][j]*self.beta[t][j]
?
return numerator/denominator
?
def cal_ksi(self, i, j, t):
"""公式 10.26"""
?
numerator = self.alpha[t][i]*self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]
denominator = 0
?
for i in range(self.N):
for j in range(self.N):
denominator += self.alpha[t][i]*self.A[i][j]*self.B[j][self.O[t+1]]*self.beta[t+1][j]
?
return numerator/denominator
?
def init(self):
"""隨機(jī)生成 A，B，Pi并保證每行相加等于 1"""
import random
for i in range(self.N):
randomlist = [random.randint(0,100) for t in range(self.N)]
Sum = sum(randomlist)
for j in range(self.N):
self.A[i][j] = randomlist[j]/Sum
?
for i in range(self.N):
randomlist = [random.randint(0,100) for t in range(self.M)]
Sum = sum(randomlist)
for j in range(self.M):
self.B[i][j] = randomlist[j]/Sum
?
def train(self, O, MaxSteps = 100):
self.T = len(O)
self.O = O
?
# 初始化
self.init()
?
step = 0
# 遞推
while step
step+=1
print(step)
tmp_A = np.zeros((self.N,self.N))
tmp_B = np.zeros((self.N,self.M))
tmp_pi = np.array([0.0]*self.N)
?
self.forward()
self.backward()
?
# a_{ij}
for i in range(self.N):
for j in range(self.N):
numerator=0.0
denominator=0.0
for t in range(self.T-1):
numerator += self.cal_ksi(i,j,t)
denominator += self.cal_gamma(i,t)
tmp_A[i][j] = numerator/denominator
?
# b_{jk}
for j in range(self.N):
for k in range(self.M):
numerator = 0.0
denominator = 0.0
for t in range(self.T):
if k == self.O[t]:
numerator += self.cal_gamma(j,t)
denominator += self.cal_gamma(j,t)
tmp_B[j][k] = numerator / denominator
?
# pi_i
for i in range(self.N):
tmp_pi[i] = self.cal_gamma(i,0)
?
self.A = tmp_A
self.B = tmp_B
self.Pi = tmp_pi
?
def generate(self, length):
import random
I = []
?
# start
ran = random.randint(0,1000)/1000.0
i = 0
while self.Pi[i]
ran -= self.Pi[i]
i += 1
I.append(i)
?
# 生成狀態(tài)序列
for i in range(1,length):
last = I[-1]
ran = random.randint(0, 1000) / 1000.0
i = 0
while self.A[last][i] < ran or self.A[last][i]<0.0001:
ran -= self.A[last][i]
i += 1
I.append(i)
?
# 生成觀測(cè)序列
Y = []
for i in range(length):
k = 0
ran = random.randint(0, 1000) / 1000.0
while self.B[I[i]][k] < ran or self.B[I[i]][k]<0.0001:
ran -= self.B[I[i]][k]
k += 1
Y.append(k)
?
return Y
?
?
?
def triangle(length):
'''三角波'''
X = [i for i in range(length)]
Y = []
?
for x in X:
x = x % 6
if x <= 3:
Y.append(x)
else:
Y.append(6-x)
return X,Y
?
?
?
def show_data(x,y):
import matplotlib.pyplot as plt
plt.plot(x, y, 'g')
plt.show()
?
return y
?
?
if __name__ == '__main__':
hmm = HMM(10,4)
tri_x, tri_y = triangle(20)
?
hmm.train(tri_y)
y = hmm.generate(100)
print(y)
x = [i for i in range(100)]
show_data(x,y)
5.2 運(yùn)行結(jié)果
5.2.1 三角波
三角波比較簡(jiǎn)單，我設(shè)置N=10，扔進(jìn)去長(zhǎng)度為20的序列，訓(xùn)練100次，下圖是其生成的長(zhǎng)度為100的序列。
可以看出效果還是很不錯(cuò)的。
5.2.2 正弦波
正弦波有些麻煩，因?yàn)橛^測(cè)序列不能太大，所以我設(shè)置N=15，M=100，扔進(jìn)去長(zhǎng)度為40的序列，訓(xùn)練100次，訓(xùn)練的非常慢，下圖是其生成的長(zhǎng)度為100的序列。
可以看出效果一般，如果改成C的代碼，并增加N應(yīng)該是能模擬的。
6 HMM的實(shí)際應(yīng)用舉例
參考文獻(xiàn)：
[2]《機(jī)器學(xué)習(xí)》周志華
[3]《統(tǒng)計(jì)學(xué)習(xí)方法》李航
延伸閱讀：
△長(zhǎng)按關(guān)注「邁微電子研發(fā)社」
△長(zhǎng)按加入「邁微電子研發(fā)社」學(xué)習(xí)輔導(dǎo)群
給個(gè)關(guān)注么么噠！
創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎(jiǎng)勵(lì)來(lái)咯，堅(jiān)持創(chuàng)作打卡瓜分現(xiàn)金大獎(jiǎng)總結(jié)
以上是生活随笔為你收集整理的python 隐马尔科夫_机器学习算法之——隐马尔可夫(Hidden Markov ModelsHMM)原理及Python实现...的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： 融e借办理条件
- 下一篇： python去年软件排行_2017年编程