python實(shí)現(xiàn)隱馬爾科夫模型HMM
更新時間：2018年03月25日 11:27:53 作者：adzhua
這篇文章主要為大家詳細(xì)介紹了python實(shí)現(xiàn)隱馬爾科夫模型HMM，具有一定的參考價值，感興趣的小伙伴們可以參考一下
一份完全按照李航<<統(tǒng)計(jì)學(xué)習(xí)方法>>介紹的HMM代碼，供大家參考，具體內(nèi)容如下
#coding=utf8
'''''
Created on 2017-8-5
里面的代碼許多地方可以精簡，但為了百分百還原公式，就沒有精簡了。
@author: adzhua
'''
import numpy as np
class HMM(object):
def __init__(self, A, B, pi):
'''''
A: 狀態(tài)轉(zhuǎn)移概率矩陣
B: 輸出觀察概率矩陣
pi: 初始化狀態(tài)向量
'''
self.A = np.array(A)
self.B = np.array(B)
self.pi = np.array(pi)
self.N = self.A.shape[0] # 總共狀態(tài)個數(shù)
self.M = self.B.shape[1] # 總共觀察值個數(shù)
# 輸出HMM的參數(shù)信息
def printHMM(self):
print ("==================================================")
print ("HMM content: N =",self.N,",M =",self.M)
for i in range(self.N):
if i==0:
print ("hmm.A ",self.A[i,:]," hmm.B ",self.B[i,:])
else:
print (" ",self.A[i,:]," ",self.B[i,:])
print ("hmm.pi",self.pi)
print ("==================================================")
# 前向算法
def forwar(self, T, O, alpha, prob):
'''''
T: 觀察序列的長度
O: 觀察序列
alpha: 運(yùn)算中用到的臨時數(shù)組
prob: 返回值所要求的概率
'''
# 初始化
for i in range(self.N):
alpha[0, i] = self.pi[i] * self.B[i, O[0]]
# 遞歸
for t in range(T-1):
for j in range(self.N):
sum = 0.0
for i in range(self.N):
sum += alpha[t, i] * self.A[i, j]
alpha[t+1, j] = sum * self.B[j, O[t+1]]
# 終止
sum = 0.0
for i in range(self.N):
sum += alpha[T-1, i]
prob[0] *= sum
# 帶修正的前向算法
def forwardWithScale(self, T, O, alpha, scale, prob):
scale[0] = 0.0
# 初始化
for i in range(self.N):
alpha[0, i] = self.pi[i] * self.B[i, O[0]]
scale[0] += alpha[0, i]
for i in range(self.N):
alpha[0, i] /= scale[0]
# 遞歸
for t in range(T-1):
scale[t+1] = 0.0
for j in range(self.N):
sum = 0.0
for i in range(self.N):
sum += alpha[t, i] * self.A[i, j]
alpha[t+1, j] = sum * self.B[j, O[t+1]]
scale[t+1] += alpha[t+1, j]
for j in range(self.N):
alpha[t+1, j] /= scale[t+1]
# 終止
for t in range(T):
prob[0] += np.log(scale[t])
def back(self, T, O, beta, prob):
'''''
T: 觀察序列的長度 len(O)
O: 觀察序列
beta: 計(jì)算時用到的臨時數(shù)組
prob: 返回值；所要求的概率
'''
# 初始化
for i in range(self.N):
beta[T-1, i] = 1.0
# 遞歸
for t in range(T-2, -1, -1): # 從T-2開始遞減；即T-2, T-3, T-4, ..., 0
for i in range(self.N):
sum = 0.0
for j in range(self.N):
sum += self.A[i, j] * self.B[j, O[t+1]] * beta[t+1, j]
beta[t, i] = sum
# 終止
sum = 0.0
for i in range(self.N):
sum += self.pi[i]*self.B[i,O[0]]*beta[0,i]
prob[0] = sum
# 帶修正的后向算法
def backwardWithScale(self, T, O, beta, scale):
'''''
T: 觀察序列的長度 len(O)
O: 觀察序列
beta: 計(jì)算時用到的臨時數(shù)組
'''
# 初始化
for i in range(self.N):
beta[T-1, i] = 1.0
# 遞歸
for t in range(T-2, -1, -1):
for i in range(self.N):
sum = 0.0
for j in range(self.N):
sum += self.A[i, j] * self.B[j, O[t+1]] * beta[t+1, j]
beta[t, i] = sum / scale[t+1]
# viterbi算法
def viterbi(self, O):
'''''
O: 觀察序列
'''
T = len(O)
# 初始化
delta = np.zeros((T, self.N), np.float)
phi = np.zeros((T, self.N), np.float)
I = np.zeros(T)
for i in range(self.N):
delta[0, i] = self.pi[i] * self.B[i, O[0]]
phi[0, i] = 0.0
# 遞歸
for t in range(1, T):
for i in range(self.N):
delta[t, i] = self.B[i, O[t]] * np.array([delta[t-1, j] * self.A[j, i] for j in range(self.N)] ).max()
phi = np.array([delta[t-1, j] * self.A[j, i] for j in range(self.N)]).argmax()
# 終止
prob = delta[T-1, :].max()
I[T-1] = delta[T-1, :].argmax()
for t in range(T-2, -1, -1):
I[t] = phi[I[t+1]]
return prob, I
# 計(jì)算gamma(計(jì)算A所需的分母；詳情見李航的統(tǒng)計(jì)學(xué)習(xí)) : 時刻t時馬爾可夫鏈處于狀態(tài)Si的概率
def computeGamma(self, T, alpha, beta, gamma):
''''''''
for t in range(T):
for i in range(self.N):
sum = 0.0
for j in range(self.N):
sum += alpha[t, j] * beta[t, j]
gamma[t, i] = (alpha[t, i] * beta[t, i]) / sum
# 計(jì)算sai(i,j)(計(jì)算A所需的分子) 為給定訓(xùn)練序列O和模型lambda時
def computeXi(self, T, O, alpha, beta, Xi):
for t in range(T-1):
sum = 0.0
for i in range(self.N):
for j in range(self.N):
Xi[t, i, j] = alpha[t, i] * self.A[i, j] * self.B[j, O[t+1]] * beta[t+1, j]
sum += Xi[t, i, j]
for i in range(self.N):
for j in range(self.N):
Xi[t, i, j] /= sum
# 輸入 L個觀察序列O，初始模型：HMM={A,B,pi,N,M}
def BaumWelch(self, L, T, O, alpha, beta, gamma):
DELTA = 0.01 ; round = 0 ; flag = 1 ; probf = [0.0]
delta = 0.0; probprev = 0.0 ; ratio = 0.0 ; deltaprev = 10e-70
xi = np.zeros((T, self.N, self.N)) # 計(jì)算A的分子
pi = np.zeros((T), np.float) # 狀態(tài)初始化概率
denominatorA = np.zeros((self.N), np.float) # 輔助計(jì)算A的分母的變量
denominatorB = np.zeros((self.N), np.float)
numeratorA = np.zeros((self.N, self.N), np.float) # 輔助計(jì)算A的分子的變量
numeratorB = np.zeros((self.N, self.M), np.float) # 針對輸出觀察概率矩陣
scale = np.zeros((T), np.float)
while True:
probf[0] =0
# E_step
for l in range(L):
self.forwardWithScale(T, O[l], alpha, scale, probf)
self.backwardWithScale(T, O[l], beta, scale)
self.computeGamma(T, alpha, beta, gamma) # (t, i)
self.computeXi(T, O[l], alpha, beta, xi) #(t, i, j)
for i in range(self.N):
pi[i] += gamma[0, i]
for t in range(T-1):
denominatorA[i] += gamma[t, i]
denominatorB[i] += gamma[t, i]
denominatorB[i] += gamma[T-1, i]
for j in range(self.N):
for t in range(T-1):
numeratorA[i, j] += xi[t, i, j]
for k in range(self.M): # M為觀察狀態(tài)取值個數(shù)
for t in range(T):
if O[l][t] == k:
numeratorB[i, k] += gamma[t, i]
# M_step。 計(jì)算pi, A, B
for i in range(self.N): # 這個for循環(huán)也可以放到for l in range(L)里面
self.pi[i] = 0.001 / self.N + 0.999 * pi[i] / L
for j in range(self.N):
self.A[i, j] = 0.001 / self.N + 0.999 * numeratorA[i, j] / denominatorA[i]
numeratorA[i, j] = 0.0
for k in range(self.M):
self.B[i, k] = 0.001 / self.N + 0.999 * numeratorB[i, k] / denominatorB[i]
numeratorB[i, k] = 0.0
#重置
pi[i] = denominatorA[i] = denominatorB[i] = 0.0
if flag == 1:
flag = 0
probprev = probf[0]
ratio = 1
continue
delta = probf[0] - probprev
ratio = delta / deltaprev
probprev = probf[0]
deltaprev = delta
round += 1
if ratio <= DELTA :
print('num iteration: ', round)
break
if __name__ == '__main__':
print ("python my HMM")
# 初始的狀態(tài)概率矩陣pi；狀態(tài)轉(zhuǎn)移矩陣A；輸出觀察概率矩陣B; 觀察序列
pi = [0.5,0.5]
A = [[0.8125,0.1875],[0.2,0.8]]
B = [[0.875,0.125],[0.25,0.75]]
O = [
[1,0,0,1,1,0,0,0,0],
[1,1,0,1,0,0,1,1,0],
[0,0,1,1,0,0,1,1,1]
]
L = len(O)
T = len(O[0]) # T等于最長序列的長度就好了
hmm = HMM(A, B, pi)
alpha = np.zeros((T,hmm.N),np.float)
beta = np.zeros((T,hmm.N),np.float)
gamma = np.zeros((T,hmm.N),np.float)
# 訓(xùn)練
hmm.BaumWelch(L,T,O,alpha,beta,gamma)
# 輸出HMM參數(shù)信息
hmm.printHMM()
以上就是本文的全部內(nèi)容，希望對大家的學(xué)習(xí)有所幫助，也希望大家多多支持腳本之家。
相關(guān)文章
這篇文章主要介紹了python如何開啟多線程問題,具有很好的參考價值,希望對大家有所幫助,如有錯誤或未考慮完全的地方,望不吝賜教2023-08-08
每個電子商務(wù)數(shù)據(jù)分析師必須掌握的一項(xiàng)數(shù)據(jù)聚類技能，如果你是一名在電子商務(wù)公司工作的數(shù)據(jù)分析師，從客戶數(shù)據(jù)中挖掘潛在價值，來提高客戶留存率很可能就是你的工作任務(wù)之一。這篇就來告訴你如何將客戶分成不同的群組，并在一段時間內(nèi)觀察每個群組的留存率2023-02-02
這篇文章主要介紹了python使用udp實(shí)現(xiàn)聊天器功能，本文通過實(shí)例代碼給大家介紹的非常詳細(xì)，具有一定的參考借鑒價值 ,需要的朋友可以參考下2018-12-12
這篇文章主要介紹了簡單的python協(xié)同過濾程序，分享了相關(guān)代碼示例，小編覺得還是挺不錯的，具有一定借鑒價值，需要的朋友可以參考下2018-01-01
python實(shí)現(xiàn)修改xml文件內(nèi)容
這篇文章主要介紹了python實(shí)現(xiàn)修改xml文件內(nèi)容，XML 指可擴(kuò)展標(biāo)記語言，是一種標(biāo)記語言，是從標(biāo)準(zhǔn)通用標(biāo)記語言（SGML）中簡化修改出來的2022-07-07
pycharm如何debug for循環(huán)里面的錯誤值(推薦)
一般debug時,在for循環(huán)里面的話,需要自己一步一步點(diǎn),如果循環(huán)幾百次那種就比較麻煩,此時可以采用try except的方式來解決,這篇文章主要介紹了pycharm如何debug for循環(huán)里面的錯誤值,需要的朋友可以參考下2024-07-07
PyQt5 關(guān)于Qt Designer的初步應(yīng)用和打包過程詳解
Qt Designer中的操作方式十分靈活，其通過拖拽的方式放置控件可以隨時查看控件效果。這篇文章主要介紹了PyQt5 關(guān)于Qt Designer的初步應(yīng)用和打包,需要的朋友可以參考下2021-09-09
如何在python中用os模塊實(shí)現(xiàn)批量移動文件
在工作中難免會遇到需要批量整理文件的情況,當(dāng)需要從一堆文件中將部分文件批量地轉(zhuǎn)移時,如果手工一一轉(zhuǎn)移難免浪費(fèi)時間,這篇文章主要給大家介紹了關(guān)于如何在python中用os模塊實(shí)現(xiàn)批量移動文件的相關(guān)資料,需要的朋友可以參考下2022-05-05
在本篇文章里小編給大家整理了關(guān)于python判斷變量是否為列表的方法，有需要的朋友們可以學(xué)習(xí)下。2020-09-09
關(guān)于python pycharm中輸出的內(nèi)容不全的解決辦法
這篇文章主要介紹了關(guān)于python pycharm中輸出的內(nèi)容不全的解決辦法,本文給大家介紹的非常詳細(xì)，具有一定的參考借鑒價值，需要的朋友可以參考下2020-01-01