首先上代碼
這里是github上的關於hmm的:鏈接
-
概率計算問題：前向-后向算法
-
學習問題：Baum-Welch算法(狀態未知)
-
預測問題：Viterbi算法
https://github.com/TimVerion/HMM_code
需要的理論基礎(可以跳過)
信息熵
首先了解一下過去化學學習的熵，熱力學中表征物質狀態的參量之一，用符號S表示，其物理意義是體系混亂程度的度量。克勞修斯於 1865 年的論文中定義了“熵” ，其中有兩句名言：“宇宙的能量是恆定的。”，“宇宙的熵趨於最大值。”
信息量：指的是一個樣本/事件所蘊含的信息，如果一個事件的概率越大，那么就可以認為該事件所蘊含的信息越少。極端情況下，比如：“太陽從東方升起”，因為是確定事件，所以不攜帶任何信息量。
信息熵：1948年，香農引入信息熵；一個系統越是有序，信息熵就越低，一個系統越是混亂，信息熵就越高，所以信息熵被認為是一個系統有序程度的度量。信息熵就是用來描述系統信息量的不確定度。
信息熵(Entropy)公式:
High Entropy(高信息熵)：表示隨機變量X是均勻分布的，各種取值情況是等概率出現的。Low Entropy(低信息熵)：表示隨機變量X各種取值不是等概率出現。可能出現有的事件概率很大，有的事件概率很小。
例子：
最大熵模型
機器學習中經常提到的最大熵的思想，就是當你要猜一個概率分布時，如果你對這個分布一無所知，那就猜熵最大的均勻分布，如果你對這個分布知道一些情況，那么，就猜滿足這些情況的熵最大的分布，就像我們使用最大似然去預測一樣。
例子：軟銀的孫正義他強調商戰要達到“不戰而屈人之兵”，避免“不打敗仗”，就得從事優勢職業。
孫正義投資雅虎，阿里，網約車火了便投資中國的滴滴和美國的Uber，也就是不把所有雞蛋放進一個籃子里，這正是最大熵原理。
與最大熵思想區分開來后，我們要知道最大熵模型就是讓信息熵最大，而所謂的條件最大熵模型，就是在一定約束下條件熵最大的模型。再直白一點就是我們要保留全部的不確定性，將風險降到最小。
也就是當從模型總體隨機抽取n組樣本觀測值后，最合理的參數估計量應該使得從模型中抽取該n組樣本觀測值的概率最大，這樣我們便構造了一個最大熵模型，沒錯這正是最大似然估計的定義。
最大熵模型，可以說是集簡與繁於一體，形式簡單，實現復雜。值得一提的是，在Google的很多產品中，比如機器翻譯，都直接或間接地用到了最大熵模型。
貝葉斯算法
三門問題：出自美國的電視游戲節目。參賽者會看見三扇關閉了的門，其中一扇的后面有一輛汽車，選中后面有車的那扇門可贏得該汽車，另外兩扇門后面則各藏有一只山羊。當參賽者選定了一扇門，但未去開啟它的時候，節目主持人開啟剩下兩扇門的其中一扇，露出其中一只山羊。主持人其后會問參賽者要不要換另一扇仍然關上的門。問題是：換另一扇門會否增加參賽者贏得汽車的機率？如果嚴格按照上述的條件，即主持人清楚地知道，自己打開的那扇門后是羊，那么答案是會。不換門的話，贏得汽車的幾率是1/3。換門的話，贏得汽車的幾率是2/3。
先驗概率P(A)：在不考慮任何情況下，A事件發生的概率
條件概率P(B|A)：A事件發生的情況下，B事件發生的概率
后驗概率P(A|B)：在B事件發生之后，對A事件發生的概率的重新評估
全概率：如果A和A'構成樣本空間的一個划分，那么事件B的概率為：A和A'的概率分別乘以B對這兩個事件的概率之和。
有上面的知識推出貝葉斯公式(后驗概率)：
貝葉斯網絡
把某個研究系統中涉及到的隨機變量，根據是否條件獨立繪制在一個有向圖中，就形成了貝葉斯網絡。貝葉斯網絡(Bayesian Network)，又稱有向無環圖模型(directed acyclicgraphical model, DAG)，是一種概率圖模型，根據概率圖的拓撲結構，考察一組隨機變量{X1,X2,...,Xn}及其N組條件概率分布(Conditional ProbabililtyDistributions, CPD)的性質
當多個特征屬性之間存在着某種相關關系的時候，使用朴素貝葉斯算法就沒法解決這類問題，那么貝葉斯網絡就是解決這類應用場景的一個非常好的算法。一般而言，貝葉斯網絡的有向無環圖中的節點表示隨機變量，可以是可觀察到的變量，或隱變量，未知參數等等。連接兩個節點之間的箭頭代表兩個隨機變量之間的因果關系(也就是這兩個隨機變量之間非條件獨立)，如果兩個節點間以一個單箭頭連接在一起，表示其中一個節點是“因”，另外一個是“果”，從而兩節點之間就會產生一個條件概率值。注意：每個節點在給定其直接前驅的時候，條件獨立於其后繼。
貝葉斯網絡的關鍵方法是圖模型，構建一個圖模型我們需要把具有因果聯系的各個變量用箭頭連在一起。貝葉斯網絡的有向無環圖中的節點表示隨機變量。連接兩個節點的箭頭代表此兩個隨機變量是具有因果關系的。貝葉斯網絡是模擬人的認知思維推理模式的，用一組條件概率以及有向無環圖對不確定性因果推理關系建模
P(a,b, c) = P(c | a,b)P(b | a)P(a)
EM算法
上面的知識點講到了MLE（最大似然估計），這里先說一下MAP(最大后驗概率)，MAP和MLE一樣，都是通過樣本估計參數θ的值；在MLE中，是使似然函數P(x|θ)最大的時候參數θ的值，MLE中假設先驗概率是一個等值的；而在MAP中，則是求θ使P(x|θ)P(θ)的值最大，這也就是要求θ值不僅僅是讓似然函數最大，同時要求θ本身出現的先驗概率也得比較大。
可以認為MAP是貝葉斯算法的一種應用:
EM算法(Expectation Maximization Algorithm, 最大期望算法)是一種迭代類型
的算法，是一種在概率模型中尋找參數最大似然估計或者最大后驗估計的算法，
其中概率模型依賴於無法觀測的隱藏變量。
EM算法流程：
-
初始化分布參數
-
重復下列兩個操作直到收斂：
E步驟：估計隱藏變量的概率分布期望函數；
M步驟：根據期望函數重新估計分布參數
EM實現的過程中還會有一個Jensen不等式知識點,它使得我們可以假設隱含數據並形成極大化模型,然后對聯合概率求最大值直到收斂。
""" 實現GMM高斯混合聚類 根據EM算法流程實現這個流程 """ import numpy as np from scipy.stats import multivariate_normal def train(x, max_iter=100): """ 進行GMM模型訓練，並返回對應的μ和σ的值(假定x數據中的簇類別數目為2) :param x: 輸入的特征矩陣x :param max_iter: 最大的迭代次數 :return: 返回一個五元組(pi, μ1， μ2，σ1，σ2) """ # 1. 獲取樣本的數量m以及特征維度n m, n = np.shape(x) # 2. 初始化相關變量 # 以每一列中的最小值作為mu1，mu1中的元素數目就是列的數目(n)個 mu1 = x.min(axis=0) mu2 = x.max(axis=0) sigma1 = np.identity(n) sigma2 = np.identity(n) pi = 0.5 # 3. 實現EM算法 for i in range(max_iter): # a. 初始化多元高斯分布（初始化兩個多元高斯混合概率密度函數） norm1 = multivariate_normal(mu1, sigma1) norm2 = multivariate_normal(mu2, sigma2) # E step # 計算所有樣本數據在norm1和norm2中的概率 tau1 = pi * norm1.pdf(x) tau2 = (1 - pi) * norm2.pdf(x) # 概率做一個歸一化操作 w = tau1 / (tau1 + tau2) # M step mu1 = np.dot(w, x) / np.sum(w) mu2 = np.dot(1 - w, x) / np.sum(1 - w) sigma1 = np.dot(w * (x - mu1).T, (x - mu1)) / np.sum(w) sigma2 = np.dot((1 - w) * (x - mu2).T, (x - mu2)) / np.sum(1 - w) pi = np.sum(w) / m # 返回最終解 return (pi, mu1, mu2, sigma1, sigma2) if __name__ == '__main__': np.random.seed(28) # 產生一個服從多元高斯分布的數據（標准正態分布的多元高斯數據） mean1 = (0, 0, 0) # x1\x2\x3的數據分布都是服從正態分布的，同時均值均為0 cov1 = np.diag((1, 1, 1)) data1 = np.random.multivariate_normal(mean=mean1, cov=cov1, size=500) # 產生一個數據分布不均衡 mean2 = (2, 2, 3) cov2 = np.array([[1, 1, 3], [1, 2, 1], [0, 0, 1]]) data2 = np.random.multivariate_normal(mean=mean2, cov=cov2, size=200) # 合並兩個數據 data = np.vstack((data1, data2)) pi, mu1, mu2, sigma1, sigma2 = train(data, 100) print("第一個類別的相關參數:") print(mu1) print(sigma1) print("第二個類別的相關參數:") print(mu2) print(sigma2) print("預測樣本屬於那個類別(概率越大就是那個類別)：") norm1 = multivariate_normal(mu1, sigma1) norm2 = multivariate_normal(mu2, sigma2) x = np.array([0, 1, 0]) print(pi * norm1.pdf(x)) # 屬於類別1的概率為:0.0275 => 0.989 print((1 - pi) * norm2.pdf(x))# 屬於類別1的概率為:0.0003 => 0.011
HMM詳解
隱馬爾可夫模型（Hidden Markov Model，HMM）作為一種統計分析模型，創立於20世紀70年代。80
年代得到了傳播和發展，成為信號處理的一個重要方向，現已成功地用於語音識別，行為識別，文字識別以及故障診斷等領域。
HMM介紹
HMM是統計模型，它用來描述一個含有隱含未知參數的馬爾可夫過程。其難點是從可觀察的參數中確定該過程的隱含參數。然后利用這些參數來作進一步的分析，例如模式識別。
馬爾可夫性質
數學定理:設{X(t), t ∈ T}是一個隨機過程，E為其狀態空間，若對於任意的t1<t2<...<tn<t,任意的x1,x2,...,xn,x∈E，隨機變量X(t)在已知變量X(t1)=x1,...,X(tn)=xn之下的條件分布函數只與X(tn)=xn有關，而與X(t1)=x1,...,X(tn-1)=xn-1無關，即條件分布函數滿足下列等式，此性質稱為馬爾可夫性；如果隨機過程滿足馬爾可夫性，則該過程稱為馬爾可夫過程。
大致意思就是我們當前狀態只受上一狀態的影響,而跟上一狀態之前的狀態無關,就叫做馬爾可夫性。
馬爾可夫鏈
馬爾可夫鏈是指具有馬爾可夫性質的隨機過程。在過程中，在給定當前信息的情況下，過去的信息狀態對於預測將來狀態是無關的。在馬爾可夫鏈的每一步，系統根據概率分布，可以從一個狀態變成另外一個狀態，也可以保持當前狀態不變。狀態的改變叫做轉移，狀態改變的相關概率叫做轉移概率。馬爾可夫鏈中的三元素是：狀態空間S、轉移概率矩陣P、初始概率分布π。
舉個例子來講：
設將天氣狀態分為晴、陰、雨三種狀態，假定某天的天氣狀態只和上一天的天氣狀態有關，狀態使用1(晴)、2(陰)、3(雨)表示，轉移概率矩陣P如下：
我們再看一下它們的各自轉移方式：
假設某天的天氣概率只與前一天的概率有關,也就是如下公式：
下面我們開始迭代，首先假設第一天的初始概率Π=[0.5,0.3,0.2]，由上式和P矩陣迭代：
這里重新假設第一天的初始概率π=[0.1,0.6,0.3]，由上式和P矩陣迭代：
由此得出規律，只要迭代次數夠多，我們初始概率π並不會影響最終概率的出現。
隱馬爾科夫模型
HMM是一種統計模型，在語音識別、行為識別、NLP、故障診斷等領域具有高效的性能。
HMM是關於時序的概率模型，描述一個含有未知參數的馬爾可夫鏈所生成的不可觀測的狀態隨機序列，再由各個狀態生成觀測隨機序列的過程。
HMM是一個雙重隨機過程---具有一定狀態的隱馬爾可夫鏈和隨機的觀測序列。HMM隨機生成的狀態隨機序列被稱為狀態序列；每個狀態成一個觀測，由此產生的觀測隨機序列，被稱為觀測序列。
HMM由隱含狀態S、可觀測狀態O、初始狀態概率矩陣π、隱含狀態轉移概率矩陣A、可觀測值轉移矩陣B(又稱為混淆矩陣，Confusion Matrix)； 
這里Π和A決定了狀態序列S，B決定了可觀測序列O，所以它們三個式整個隱馬的關鍵：
這里S和O是相對於整個集合而言的，我們一般只會度量前T個狀態序列I和觀測序列Q:
這里每個q發生的概率只受i的影響,而每一個i只受前一個i的影響。
HMM案例
假設有三個盒子，編號為1,2,3；每個盒子都裝有黑白兩種顏色的小球，球的比例如下：
| 編號 | 白球 | 黑球 |
|---|---|---|
| 1 | 4 | 6 |
| 2 | 8 | 2 |
| 3 | 5 | 5 |
按照下列規則的方式進行有放回的抽取小球，得到球顏色的觀測序列：
-
按照π的概率選擇一個盒子，從盒子中隨機抽取出一個小球(B)，記錄顏色后，放回盒子中；
-
按照某種條件概率(A)選擇新的盒子，重復該操作；
-
最終得到觀測序列：“白黑白白黑“
狀態集合：S={盒子1，盒子2，盒子3}觀測集合：O={白，黑}設狀態序列和觀測序列的長度T=5，並假設如下A，B，π
狀態轉移概率矩陣A:
觀測概率矩陣B:
初始概率分布π:
這時候我們面對自己假設的A，B，π，那它們得到觀測序列“白黑白白黑“的概率是多少？
現在假設條件下，唯一未知得就是狀態序列，所以這里就要嘗試使用不同的狀態序列，並找到在一個在估計模型下的擁有最大概率的狀態序列。
有了狀態序列我們就可以預測再次抽取時發生事件的概率，在這個例子里就是我們每次到底抽的時那個盒子，但是估計模型是假設的所以在去求狀態序列的前面，應該先去求得最優得A，B，π使得在這個模型下觀測序列發生得概率最大。
這里便出現了兩個問題：
1，如何求得一個估計模型，讓觀測序列在該模型下發生概率最大。
2，如何求得隱藏得狀態序列，讓它在最優估計模型下發生得概率最大。
這里為什么要求概率最大，便涉及到了最大熵模型。
應用到實際便成了三個問題，就是多一個我們如何該觀測序列在估計模型上出現得概率。
三個問題
概率計算問題：前向-后向算法
給定模型λ=(A,B,π)和觀測序列Q={q1,q2,...,qT}，計算模型λ下觀測到序列Q出現的概率P(Q|λ)
學習問題：Baum-Welch算法(狀態未知)
已知觀測序列Q={q1,q2,...,qT}，估計模型λ=(A,B,π)的參數，使得在該模型下觀測序列P(Q|λ)最大。
預測問題：Viterbi算法
給定模型λ=(A,B,π)和觀測序列Q={q1,q2,...,qT}，求給定觀測序列條件概率P(I|Q，λ)最大的狀態序列I
概率計算問題（前后向）
前向概率-后向概率指的其實是在一個觀測序列中，時刻t對應的狀態為si的概率值轉換過來的信息。
所以這兩種都可以解決概率計算問題，也就是看懂一個就可以繼續下一個問題。
首先寫出在狀態序列為s_i的情況下，觀測序列出現的概率。
前向算法定義：給定λ（A,B,π），定義到時刻t部分觀測序列為q1,q2,...,qt且狀態為si的概率為前向概率。此時我們加設β為1,所以表達方程記做：
-
初值
先計算第一個α：
-
遞推
求出第一個后,遞推使t = 1,2,...,T-1
-
結果
這樣狀態序列為s_i的情況下，觀測序列出現的概率就變成了
-
[前向傳播代碼] https://github.com/TimVerion/HMM_code/blob/master/hmm/forward_probability.py 前向傳播代碼
-
# 偽代碼，不可運行 # 更新初值(t=1) for i in n_range: alpha[0][i] = pi[i] * B[i][fetch_index_by_obs_seq_f(Q, 0)] # 迭代更新其它時刻 T = len(Q) tmp = [0 for i in n_range] for t in range(1, T): for i in n_range: # 1. 計算上一個時刻t-1累積過來的概率值 for j in n_range: tmp[j] = alpha[t - 1][j] * A[j][i] # 2. 更新alpha的值 alpha[t][i] = np.sum(tmp) * B[i][fetch_index_by_obs_seq_f(Q, t)]
后向傳播定義:給定λ，定義到時刻t狀態為si的前提下，從t+1到T部分觀測序列為qt+1,qt+2,...,qT的概率為后向概率。記做:
1. 初值
先計算第一個β,在我們進行前向算法的時候我們假設β為1：
-
遞推
求出第一個后,遞推使t = T-1,T-2...,1
-
結果
這樣狀態序列為s_i的情況下，觀測序列出現的概率就變成了:
[ 后向傳播代碼] https://github.com/TimVerion/HMM_code/blob/master/hmm/backward_probability.py 后向傳播代碼
# 更新初值(t=T) for i in n_range: beta[T - 1][i] = 1 # 迭代更新其它時刻 tmp = [0 for i in n_range] for t in range(T - 2, -1, -1): for i in n_range: # 1. 計算到下一個時刻t+1的概率值 for j in n_range: tmp[j] = A[i][j] * beta[t + 1][j] * B[j][fetch_index_by_obs_seq_f(Q, t +1)] # 2. 更新beta的值 beta[t][i] = np.sum(tmp)
學習問題(鮑姆韋爾奇)
根據概率計算問題我們可以求的兩個可以幫助我們解決學習問題的概率:
1,單個狀態的概率
將給定模型λ和觀測序列Q的情況下，在時刻t處於狀態si的概率，記做：
單個狀態概率的意義主要是用於判斷在每個時刻最可能存在的狀態，從而可以得到一個狀態序列作為最終的預測結果。
推導過程：
2,兩個狀態的聯合概率
將給定模型λ和觀測序列Q的情況下，在時刻t處於狀態si並且在t+1時刻處於狀態sj的概率，記做：
也就是：
3,解決學習問題
若訓練數據包含觀測序列和狀態序列，則HMM的學習問題非常簡單，是監督學習算法。我們甚至在知道觀測序列和狀態序列的時候我們可以利用大數定理求出最優模型。
但是我們的訓練數據只包含觀測序列，則HMM的學習問題需要使用EM算法求解，是非監督學習算法,這里使用的EM算法也叫做鮑姆韋爾奇。
那么我們想一下，這里假設所有的觀測數據為Q={q1,q2,...,qT}，所有的隱狀態為I={i1,i2,...,iT}，則完整的數據為(O,I)，完整數據的對數似然函數為ln(p(Q,I;λ)); 然后直接使用EM算法的方式就可以來進行參數估計了。（EM算法不懂見上）
-
初始化分布參數
-
重復下列兩個操作直到收斂：
E步驟：估計隱藏變量的概率分布期望函數:
M步驟：根據期望函數重新估計分布參數
π的求解
當我們了解L函數,也就是我們需要極大化的函數,可以直接對π求導,然后讓偏導等於零。
這里寫到γ是不是很熟悉,沒錯就是我們是上面求的單個狀態的概率。
A的求解
和上一個變量求解一樣，我們需要極大化L
在上面的式子里，我們求aij也就是A的時候，我們必須使用單個狀態的概率和兩個狀態的聯合概率。
B的求解
同上我們對B求偏導，然后讓偏導等於零：
經過極大化L函數我們可以求得π、a、b的值,這樣我們便得到了一個最優模型。
# 1. 迭代更新(EM算法思想類型) for time in range(max_iter): # a. 在當前的pi，A，B的情況下對觀測序列Q分別計算alpha、beta、gamma和ksi forward.calc_alpha(pi, A, B, Q, alpha, fetch_index_by_obs_seq_f) backward.calc_beta(pi, A, B, Q, beta, fetch_index_by_obs_seq_f) single.calc_gamma(alpha, beta, gamma) continuous.calc_ksi(alpha, beta, A, B, Q, ksi, fetch_index_by_obs_seq # b. 更新pi、A、B的值 # b.1. 更新pi值 for i in n_range: pi[i] = gamma[0] # b.2. 更新狀態轉移矩陣A的值 tmp1 = np.zeros(T - 1) tmp2 = np.zeros(T - 1) for i in n_range: for j in n_range: # 獲取所有時刻從狀態i轉移到狀態j的值 for t in t_1_range: tmp1[t] = ksi[t][i][j] tmp2[t] = gamma[t] # 更新狀態i到狀態j的轉移概率 A[i][j] = np.sum(tmp1) / np.sum(tm # b.3. 更新狀態和觀測值之間的轉移矩陣 for i in n_range: for k in m_range: tmp1 = np.zeros(T) tmp2 = np.zeros(T) # 獲取所有時刻從狀態i轉移到觀測值k的概率和 number = 0 for t in t_range: if k == fetch_index_by_obs_seq_f(Q, t): 如果序列Q中時刻t對應的觀測值就是k，那么進行統計這個時刻t為狀態i的概率值 tmp1[t] = gamma[t][i] number + tmp2[t] = gamma[t] # 更新狀態i到觀測值k之間的轉移概率 if number == 0: # 沒有轉移，所以為0 B[i][k] = 0 else: # 有具體值，那么進行更新操作 B[i][k] = np.sum(tmp1) / np.sum(tem2)
預測問題（維特比）
Viterbi算法實際是用動態規划的思路求解HMM預測問題，求出概率最大的“路徑”，每條“路徑”對應一個狀態序列。
動態規划是運籌學的一個分支，是求解決策過程最優化的數學方法。把多階段過程轉化為一系列單階段問題，利用各階段之間的關系，逐個求解。
比如斐波那契數列和背包問題都利用了這個思想。
# 1. 計算t=1的時候delta的值 for i in n_range: delta[0][i] = pi[i] * B[i][fetch_index_by_obs_seq_f(Q, 0)] # 2. 更新其它時刻的值 for t in range(1, T): for i in n_range: # 當前時刻t的狀態為i # a. 獲取最大值 max_delta = -1 for j in n_range: # j表示的是上一個時刻的狀態值 tmp = delta[t - 1][j] * A[j][i] if tmp > max_delta: max_delta = tmp pre_index[t][i] = j # b. 更新值 delta[t][i] = max_delta * B[i][fetch_index_by_obs_seq_f(Q, t)] # 3. 解碼操作，查找到最大的結果值 decode = [-1 for i in range(T)]
隱馬爾可夫模型(HMM)是由馬爾可夫過程衍生出的概率圖模型，常被用於語音模式識別、生物基因序列標記、金融時間序列預測等。
HMM在我們生活中無處不在，舉個簡單的例子：
身邊的朋友一般出行有騎共享單車出行的，說下雨就比晴天騎車的人數少，當我們聽朋友說今天外面共享單車被騎沒了，我們可以推斷出今天天氣不錯，這里顯式狀態是出行，而隱狀態是天氣。
所以我們很早就已經會HMM算法的思想了，只是差這幾個算法和證明。
這里說的HMM不是買零食的韓梅梅
馬爾可夫模型解決語音識別
時間點 t的隱藏條件和時間點 t-1的隱藏條件有關。
因為人類語音擁有前后的關系，可以從語義與發音兩點來看：
單字的發音擁有前后關系：例如"They are"常常發音成"They're"，或是"Did you"會因為"you"的發音受"did"的影響，常常發音成"did ju"，而且語音識別中用句子的發音來進行分析，因此需要考慮到每個音節的前后關系，才能夠有較高的准確率。
句子中的單字有前后關系：從英文文法來看，主詞后面常常接助動詞或是動詞，動詞后面接的會是受詞或介系詞。而或是從單一單字的使用方法來看，對應的動詞會有固定使用的介系詞或對應名詞。因此分析語音頻息時需要為了提升每個單字的准確率，也需要分析前后的單字。
馬爾可夫模型將輸入消息視為一單位一單位，接着進行分析，與人類語音模型的特性相似。語音系統識別的單位為一個單位時間內的聲音。利用梅爾倒頻譜等語音處理方法，轉換成一個發音單位，為離散型的信息。而馬爾可夫模型使用的隱藏條件也是一個個被數據包的 x(t)，因此使用馬爾可夫模型來處理聲音頻號比較合適。