Python實現HMM（隱馬爾可夫模型）
前幾天用MATLAB實現了HMM的程式碼，這次用python寫了一遍，依據仍然是李航博士的《統計學習方法》
由於第一次用python，所以程式碼可能會有許多缺陷，但是所有程式碼都用書中的例題進行了測試，結果正確。
這裡想說一下python，在編寫HMM過程中參看了之前寫的MATLAB程式，發現他們有太多相似的地方，用到了numpy庫，在python程式碼過程中最讓我頭疼的是陣列角標，和MATLAB矩陣角標從1開始不同，numpy庫陣列角標都是從0開始，而且陣列的維數也需要謹慎，一不小心就會出現too many indices for array的錯誤。程式中最後是維特比演算法，在執行過程中出現了__main__:190: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future的警告，還沒有去掉這個警告，查了一下說不影響結果，後面會去解決這個問題，下面貼出我的程式碼
# -*- coding: utf-8 -*- """ Created on Thu Feb 16 19:28:39 2017 2017-4-2 ForwardBackwardAlg函式功能：實現前向演算法 理論依據：李航《統計學習方法》 2017-4-5 修改了ForwardBackwardAlg函式名稱為ForwardAlgo以及輸出的alpha陣列形式 完成了BackwardAlgo函式功能：後向演算法 以及函式FBAlgoAppli：計算在觀測序列和模型引數確定的情況下， 某一個隱含狀態對應相應的觀測狀態的概率 2017-4-6 完成BaumWelchAlgo函式一次迭代 2017-4-7 實現維特比演算法 @author: sgp """ import numpy as np #輸入格式如下： #A = np.array([[.5,.2,.3],[.3,.5,.2],[.2,.3,.5]]) #B = np.array([[.5,.5],[.4,.6],[.7,.3]]) #Pi = np.array([[.2,.4,.4]]) #O = np.array([[1,2,1]]) #應用ndarray在陣列之間進行相互運算時，一定要確保陣列維數相同！ #比如： #In[93]:m = np.array([1,2,3,4]) #In[94]:m #Out[94]: array([1, 2, 3, 4]) #In[95]:m.shape #Out[95]: (4,) #這裡表示的是一維陣列 #In[96]:m = np.array([[1,2,3,4]]) #In[97]:m #Out[97]: array([[1, 2, 3, 4]]) #In[98]:m.shape #Out[98]: (1, 4) #而這裡表示的就是二維陣列 #注意In[93]和In[96]的區別,多一對中括號！！ #N = A.shape[0]為陣列A的行數， H = O.shape[1]為陣列O的列數 #在下列各函式中，alpha陣列和beta陣列均為N*H二維陣列，也就是橫向座標是時間，縱向是狀態 def ForwardAlgo(A,B,Pi,O): N = A.shape[0]#陣列A的行數 M = A.shape[1]#陣列A的列數 H = O.shape[1]#陣列O的列數 sum_alpha_1 = np.zeros((M,N)) alpha = np.zeros((N,H)) r = np.zeros((1,N)) alpha_1 = np.multiply(Pi[0,:], B[:,O[0,0]-1]) alpha[:,0] = np.array(alpha_1).reshape(1,N)#alpha_1是一維陣列，在使用np.multiply的時候需要升級到二維陣列。#錯誤是IndexError: too many indices for array
其實程式碼就是翻譯的公式，李航博士的書中有比較詳細的推理過程，或者去找一些專業的論文文獻進一步瞭解，這裡僅僅是實現了最簡單的應用，其應用例項如下：for h in range(1,H): for i in range(N): for j in range(M): sum_alpha_1[i,j] = alpha[j,h-1] * A[j,i] r = sum_alpha_1.sum(1).reshape(1,N)#同理，將陣列升級為二維陣列 alpha[i,h] = r[0,i] * B[i,O[0,h]-1] #print("alpha矩陣: \n %r" % alpha) p = alpha.sum(0).reshape(1,H) P = p[0,H-1] #print("觀測概率: \n %r" % P) #return alpha return alpha, P def BackwardAlgo(A,B,Pi,O): N = A.shape[0]#陣列A的行數 M = A.shape[1]#陣列A的列數 H = O.shape[1]#陣列O的列數 #beta = np.zeros((N,H)) sum_beta = np.zeros((1,N)) beta = np.zeros((N,H)) beta[:,H-1] = 1 p_beta = np.zeros((1,N)) for h in range(H-1,0,-1): for i in range(N): for j in range(M): sum_beta[0,j] = A[i,j] * B[j,O[0,h]-1] * beta[j,h] beta[i,h-1] = sum_beta.sum(1) #print("beta矩陣: \n %r" % beta) for i in range(N): p_beta[0,i] = Pi[0,i] * B[i,O[0,0]-1] * beta[i,0] p = p_beta.sum(1).reshape(1,1) #print("觀測概率: \n %r" % p[0,0]) return beta, p[0,0] def FBAlgoAppli(A,B,Pi,O,I): #計算在觀測序列和模型引數確定的情況下，某一個隱含狀態對應相應的觀測狀態的概率 #例題參考李航《統計學習方法》P189習題10.2 #輸入格式： #I為二維陣列，存放所求概率P(it = qi,O|lambda)中it和qi的角標t和i，即P=[t,i] alpha,p1 = ForwardAlgo(A,B,Pi,O) beta,p2 = BackwardAlgo(A,B,Pi,O) p = alpha[I[0,1]-1,I[0,0]-1] * beta[I[0,1]-1,I[0,0]-1] / p1 return p def GetGamma(A,B,Pi,O): N = A.shape[0]#陣列A的行數 H = O.shape[1]#陣列O的列數 Gamma = np.zeros((N,H)) alpha,p1 = ForwardAlgo(A,B,Pi,O) beta,p2 = BackwardAlgo(A,B,Pi,O) for h in range(H): for i in range(N): Gamma[i,h] = alpha[i,h] * beta[i,h] / p1 return Gamma def GetXi(A,B,Pi,O): N = A.shape[0]#陣列A的行數 M = A.shape[1]#陣列A的列數 H = O.shape[1]#陣列O的列數 Xi = np.zeros((H-1,N,M)) alpha,p1 = ForwardAlgo(A,B,Pi,O) beta,p2 = BackwardAlgo(A,B,Pi,O) for h in range(H-1): for i in range(N): for j in range(M): Xi[h,i,j] = alpha[i,h] * A[i,j] * B[j,O[0,h+1]-1] * beta[j,h+1] / p1 #print("Xi矩陣: \n %r" % Xi) return Xi def BaumWelchAlgo(A,B,Pi,O): N = A.shape[0]#陣列A的行數 M = A.shape[1]#陣列A的列數 Y = B.shape[1]#陣列B的列數 H = O.shape[1]#陣列O的列數 c = 0 Gamma = GetGamma(A,B,Pi,O) Xi = GetXi(A,B,Pi,O) Xi_1 = Xi.sum(0) a = np.zeros((N,M)) b = np.zeros((M,Y)) pi = np.zeros((1,N)) a_1 = np.subtract(Gamma.sum(1),Gamma[:,H-1]).reshape(1,N) for i in range(N): for j in range(M): a[i,j] = Xi_1[i,j] / a_1[0,i] #print(a) for y in range(Y): for j in range(M): for h in range(H): if O[0,h]-1 == y: c = c + Gamma[j,h] gamma = Gamma.sum(1).reshape(1,N) b[j,y] = c / gamma[0,j] c = 0 #print(b) for i in range(N): pi[0,i] = Gamma[i,0] #print(pi) return a,b,pi def BaumWelchAlgo_n(A,B,Pi,O,n):#計算迭代次數為n的BaumWelch演算法 for i in range(n): A,B,Pi = BaumWelchAlgo(A,B,Pi,O) return A,B,Pi def viterbi(A,B,Pi,O): N = A.shape[0]#陣列A的行數 M = A.shape[1]#陣列A的列數 H = O.shape[1]#陣列O的列數 Delta = np.zeros((M,H)) Psi = np.zeros((M,H)) Delta_1 = np.zeros((N,1)) I = np.zeros((1,H)) for i in range(N): Delta[i,0] = Pi[0,i] * B[i,O[0,0]-1] for h in range(1,H): for j in range(M): for i in range(N): Delta_1[i,0] = Delta[i,h-1] * A[i,j] * B[j,O[0,h]-1] Delta[j,h] = np.amax(Delta_1) Psi[j,h] = np.argmax(Delta_1)+1 print("Delta矩陣: \n %r" % Delta) print("Psi矩陣: \n %r" % Psi) P_best = np.amax(Delta[:,H-1]) psi = np.argmax(Delta[:,H-1]) I[0,H-1] = psi + 1 for h in range(H-1,0,-1): I[0,h-1] = Psi[I[0,h]-1,h] print("最優路徑概率: \n %r" % P_best) print("最優路徑: \n %r" % I)
輸入資料格式：
In[117]:A
Out[117]:
array([[ 0.5, 0.2, 0.3],
[ 0.3, 0.5, 0.2],
[ 0.2, 0.3, 0.5]])
In[118]:B
Out[118]:
array([[ 0.5, 0.5],
[ 0.4, 0.6],
[ 0.7, 0.3]])
In[119]:Pi
Out[119]: array([[ 0.2, 0.4, 0.4]])
In[120]:O
Out[120]: array([[1, 2, 1]])
輸出結果為：
In[101]:alpha,p = ForwardAlgo(A,B,Pi,O)
In[102]:alpha
Out[102]:
array([[ 0.1 , 0.077 , 0.04187 ],
[ 0.16 , 0.1104 , 0.035512],
[ 0.28 , 0.0606 , 0.052836]])
In[103]:p
Out[103]: 0.130218
In[104]:beta,p1 = BackwardAlgo(A,B,Pi,O)
In[105]:beta
Out[105]:
array([[ 0.2451, 0.54 , 1. ],
[ 0.2622, 0.49 , 1. ],
[ 0.2277, 0.57 , 1. ]])
In[106]:p1
Out[106]: 0.130218
In[107]:gamma = GetGamma(A,B,Pi,O)
In[108]:gamma
Out[108]:
array([[ 0.18822283, 0.31931069, 0.32153773],
[ 0.32216744, 0.41542644, 0.27271191],
[ 0.48960973, 0.26526287, 0.40575036]])
In[109]:xi = GetXi(A,B,Pi,O)
In[110]:xi
Out[110]:
array([[[ 0.1036723 , 0.04515505, 0.03939548],
[ 0.09952541, 0.18062019, 0.04202184],
[ 0.11611298, 0.1896512 , 0.18384555]],
[[ 0.14782903, 0.04730529, 0.12417638],
[ 0.12717136, 0.16956181, 0.11869327],
[ 0.04653735, 0.05584481, 0.16288071]]])
In[111]:a,b,pi = BaumWelchAlgo_n(A,B,Pi,O,5)
In[112]:a
Out[112]:
array([[ 0.43972438, 0.15395857, 0.40631705],
[ 0.309058 , 0.45055446, 0.24038754],
[ 0.3757005 , 0.50361975, 0.12067975]])
In[113]:b
Out[113]:
array([[ 0.50277235, 0.49722765],
[ 0.49524289, 0.50475711],
[ 0.91925551, 0.08074449]])
In[114]:pi
Out[114]: array([[ 0.08435301, 0.18040718, 0.73523981]])
In[115]:viterbi(A,B,Pi,O)
Delta矩陣:
array([[ 0.1 , 0.028 , 0.00756],
[ 0.16 , 0.0504 , 0.01008],
[ 0.28 , 0.042 , 0.0147 ]])
Psi矩陣:
array([[ 0., 3., 2.],
[ 0., 3., 2.],
[ 0., 3., 3.]])
最優路徑概率:
0.014699999999999998
最優路徑:
array([[ 3., 3., 3.]])
__main__:192: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
相關推薦
Python實現HMM（隱馬爾可夫模型）
前幾天用MATLAB實現了HMM的程式碼，這次用python寫了一遍，依據仍然是李航博士的《統計學習方法》 由於第一次用python，所以程式碼可能會有許多缺陷，但是所有程式碼都用書中的例題進行了測試，結果正確。 這裡想說一下python，在編寫HMM過程中參看了之前寫的M
快速通俗理解HMM（隱馬爾可夫模型）
導讀 本文的結構安排如下目錄。首先介紹了HMM的基本定義；然後引出HMM的3個基本問題中；在前兩個環節中通過一個淺顯易懂的案例幫助理解；其次是針對3個問題中用到演算法的理論介紹；最後介紹了預測天氣的HMM經典案例，並附有python程式碼。 目錄 案
一文搞懂HMM（隱馬爾可夫模型）
什麼是熵(Entropy) 簡單來說，熵是表示物質系統狀態的一種度量，用它老表徵系統的無序程度。熵越大，系統越無序，意味著系統結構和運動的不確定和無規則；反之，，熵越小，系統越有序，意味著具有確定和有規則的運動狀態。熵的中文意思是熱量被溫度除的商。負熵是物質系統有序化，組織化，複雜化狀態的一種度量。 熵最早來
HMM（隱馬爾科夫模型）
一、定義 從定義可以看出，隱馬爾可夫模型做了兩個基本假設： 二、 隱馬爾可夫模型的三要素：初始狀態概率矩陣、狀態轉移概率矩陣A、觀測概率矩陣 B 三、隱馬爾科夫的3個基本問題 3.1 概率計算演算法 3.1.1 直
基於HMM（隱馬爾科夫鏈）和DNN的歌聲合成系統成果樣本
基於HMM（隱馬爾科夫鏈）和DNN的歌聲合成系統成果樣本 閒話少說 樣本連線 閒話少說 最近工作的一點成績。 樣本連線 日文歌曲1. 日文歌曲2.； 日文歌曲3.； 日文歌曲4.；
隱馬爾可夫模型（《統計學習方法》、python實現）
本文是《統計學習方法》第10章的筆記，用一段167行的Python程式碼實現了隱馬模型觀測序列的生成、前向後向演算法、Baum-Welch無監督訓練、維特比演算法。公式與程式碼相互對照，循序漸進。HMM算是個特別常見的模型，早在我沒有挖ML這個坑的時候，就已經在用HMM做基於
ML--HMM(隱馬爾可夫模型及python的實現2)
1.HMM的應用1，這個程式碼不知道出處了，若有侵權請聯絡本文作者刪除，註釋為本人所加。 2.對基本的HMM需要進一步瞭解的，請戳這裡 3.下面是HMM程式碼的解釋之一 # _*_ codin
隱馬爾可夫模型（HMM）和 jieba分詞原始碼的理解
在理解隱馬爾可夫模型（HMM）時，看到的很好的部落格，記錄一下： 1. 隱馬爾可夫模型(HMM) - 1 - 基本概念：http://blog.csdn.net/xueyingxue001/article/details/51435728 2.隱馬爾可夫模型(HMM) - 2 -
隱馬爾可夫模型（HMM）及Viterbi演算法
HMM簡介 對於演算法愛好者來說，隱馬爾可夫模型的大名那是如雷貫耳。那麼，這個模型到底長什麼樣？具體的原理又是什麼呢？有什麼具體的應用場景呢？本文將會解答這些疑惑。 本文將通過具體形象的例子來引入該模型，並深入探究隱馬爾可夫模型及Viterbi演算法，希望能對大家有所啟發。
ml課程：概率圖模型—貝葉斯網路、隱馬爾可夫模型相關（含程式碼實現）
以下是我的學習筆記，以及總結，如有錯誤之處請不吝賜教。 本文主要介紹機器學習中的一個分支——概率圖模型、相關基礎概念以及樸素貝葉斯、隱馬爾可夫演算法，最後還有相關程式碼案例。 說到機器學習的起源，可以分為以下幾個派別： 連線主義：又稱為仿生學派(bionicsism)或生理學派
HMM隱馬爾可夫模型（HMM）攻略
隱馬爾可夫模型 (Hidden Markov Model，HMM) 最初由 L. E. Baum 和其它一些學者發表在一系列的統計學論文中，隨後在語言識別，自然語言處理以及生物資訊等領域體現了很大的價值。平時，經常能接觸到涉及 HMM 的相關文章，一直沒有仔細研究過，都
隱馬爾可夫模型（Hidden Markov Model, HMM)
1. 基本概念 1.1 HMM定義 隱馬爾可夫模型（簡稱HMM），描述由隱藏的馬爾可夫鏈隨機生成不可觀測的狀態序列，再由各個狀態生成觀測序列的過程。隱藏的馬爾可夫鏈生成的不可觀測的狀態隨機序列稱為狀態序列。每一個狀態可生成一個觀測，各個狀態產生的隨機序列稱為觀測
統計學習方法_隱馬爾可夫模型HMM實現
這裡用到的資料集是三角波，使用長度20的序列訓練100次，生成長度為100的序列。HMM的初始化非常重要，這裡採用隨機初始化。 #!/usr/bin/env python3 # -*- coding: utf-8 -*- import csv import random
隱馬爾可夫模型（HMM）和Viterbi演算法
1. 隱馬爾可夫模型（HMM） 在說隱馬爾可夫模型前還有一個概念叫做“馬爾科夫鏈”，既是在給定當前知識或資訊的情況下，觀察物件過去的歷史狀態對於預測將來是無關的。也可以說在觀察一個系統變化的時候，他的下一個狀態如何的概率只需要觀察和統計當前的狀態即可正確得出。隱馬爾可夫鏈和貝葉
隱馬爾可夫模型（HMM）詳解
上邊的圖示都強調了 HMM 的狀態變遷。而下圖則明確的表示出模型的演化，其中綠色的圓圈表示隱藏狀態，紫色圓圈表示可觀察到狀態，箭頭表示狀態之間的依存概率，一個 HMM 可用一個5元組 { N, M, π，A，B } 表示，其中 N 表示隱藏狀態的數量，我們要麼知道確切的值，要麼猜測該值，M 表示可觀測狀態的數
隱馬爾可夫模型HMM（二）概率計算問題
摘自 1.李航的《統計學習方法》 2.http://www.cnblogs.com/pinard/p/6955871.html 一、概率計算問題 上一篇介紹了概率計算問題是給定了λ（A，B，π），計算一個觀測序列O出現的概率，即求P（O|λ）。 用三種方法，直接計演算法，前向演算法，
隱馬爾可夫模型（HMM）攻略
上邊的圖示都強調了 HMM 的狀態變遷。而下圖則明確的表示出模型的演化，其中綠色的圓圈表示隱藏狀態，紫色圓圈表示可觀察到狀態，箭頭表示狀態之間的依存概率，一個 HMM 可用一個5元組 { N, M, π，A，B } 表示，其中 N 表示隱藏狀態的數量，我們要麼知道確切的值，要麼猜測該值，M 表示可觀測狀態的
隱馬爾可夫模型（HMM）
隱馬爾可夫模型原理部分可以概括為三句話：一個定義、兩個假設、三個問題 HMM是一個五元組（Y，X，π，A，B），其中Y是狀態（輸出）的集合，X是觀察值（輸入）集合，π是初始狀態的概率，A是狀態
HMM隱馬爾可夫模型來龍去脈（一）
目錄 隱馬爾可夫模型HMM學習導航 一、認識貝葉斯網路 1、概念原理介紹 2、舉例解析 二、馬爾可夫模型 1、概念原理介紹 2、舉例解析 三、隱馬爾可夫模型 1、概念原理介紹 2、舉例解析 四、隱馬爾可夫模型簡單實現 五、完整程式碼 六、
HMM隱馬爾可夫模型來龍去脈（二）
目錄 前言 預備知識 一、估計問題 1、問題推導 2、前向演算法/後向演算法 二、序列問題 &