hmm隱馬爾可夫真的那麼難嗎？

時間 2019-12-11

標籤 hmm 那麼简体版

原文原文鏈接

hmm隱馬爾可夫真的那麼難嗎？

首先上代碼

這裏是github上的關於hmm的:連接git

機率計算問題：前向-後向算法github
學習問題：Baum-Welch算法(狀態未知)算法
預測問題：Viterbi算法網絡

https://github.com/TimVerion/HMM_codedom

須要的理論基礎(能夠跳過)

信息熵

首先了解一下過去化學學習的熵，熱力學中表徵物質狀態的參量之一，用符號S表示，其物理意義是體系混亂程度的度量。克勞修斯於 1865 年的論文中定義了「熵」，其中有兩句名言：「宇宙的能量是恆定的。」，「宇宙的熵趨於最大值。」機器學習

信息量：指的是一個樣本/事件所蘊含的信息，若是一個事件的機率越大，那麼就能夠認爲該事件所蘊含的信息越少。極端狀況下，好比：「太陽從東方升起」，由於是肯定事件，因此不攜帶任何信息量。ide

信息熵：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)]