Hidden Markov Model

#HMM隱馬爾科夫模型 ###①通俗的理解 首先舉一個例子，扔骰子，有三種骰子，第一個是比較常見的6個面，每個面的機率都是1/6。第二個只有4個面，，每個面的機率是1/4。第三個有8個面，，每個面的機率是1/8。

首先先選擇一個骰子，挑到每個骰子的機率1/3，而後投擲，能夠獲得 。最後會獲得一堆的序列，好比會獲得 等等，這種序列叫可見狀態序列，但在HMM裏面，還存在一個隱含狀態鏈，好比這個狀態鏈多是 。從8個面出來的，6個面投出來的。因此隱馬爾科夫模型裏面有兩串序列，可觀測序列和隱含狀態序列。 通常來說，說到隱馬爾科夫鏈其實就是隱含的狀態序列了，markov鏈各個元素直接是存在了轉化關係的，好比一開始是D6下一個狀態是D4，D6，D8的機率都是1/3，這樣設置平均值只是爲了初始化而已，通常這種值實際上是能夠隨意改變的，好比這裏下一個狀態序列是D8，那麼就能夠設置轉化到D6的機率是0.1，D4的機率是0.1，D8的機率是0.8。這樣就又是一個很新的HMM了。 一樣的，儘管可見的狀態之間沒有轉換機率，可是隱含狀態和可見狀態之間有一個機率叫作輸出機率，也叫發射機率，就這個例子來講，D6輸出1的機率就是6了，產生23456的機率都是1/6。 因此，HMM首先是有兩個序列，可觀測序列，隱含序列。每個隱含序列裏面的元素直接是存在着轉化機率的，也就是用機率來描述轉化到下一個狀態的可能性；而隱含序列和觀測序列之間也有一個機率，發射機率，也叫輸出機率。

隱含狀態關係圖：

每個箭頭都會對應着一個機率，轉換機率。 若是知道了轉換機率和隱含機率那麼求解觀測序列會很簡單，但若是是這樣那麼這個模型就沒有什麼能夠研究的了，因此確定是會缺失一點值的，好比知道觀測序列和參數，可是不知道隱含序列是什麼，或者是知道了 觀測序列想求參數，又或者是知道了這個參數和隱含序列求這個觀測序列出現的參數是什麼。其實這些就對應了後面的三個問題。HMM的核心其實也就是這三個問題。 ###②三個基本問題通俗解釋 1.知道骰子有幾種（知道隱含序列有多少種），也就是幾種骰子，上面就有三種，每種骰子是什麼（轉換機率），根據骰子扔出的結果（可觀測狀態鏈），想知道每次投的是哪一種骰子（隱含狀態序列）？ 這問題其實就是解碼問題了，知道觀測序列求隱含序列。其實有兩種解法，一種就是暴力求解：其實就是maximum likelihood所有乘起來一塊兒算便可。另外一種就厲害帶了，遞推的方法，前向算法，後向算法，前向-後向算法。 2.仍是知道骰子有幾種隱含狀態（隱含狀態的數量），知道骰子是什麼（知道轉換機率），根據骰子投出的結果我想知道投出這個結果的機率是多少？ 這個問題看起來貌似沒有什麼太大的做用，可是其實是用於對於模型的檢測，若是咱們投出來的結果都對應了很小的機率，那麼模型可能就是錯誤的，有人換了咱們的骰子。 3.知道骰子有幾種（知道隱含序列的種類），不知道每種骰子是什麼（不知道轉換機率），實驗獲得屢次骰子的結果（可觀測序列），如今想知道每種骰子是什麼（轉換機率） 這很重要，後面會用到來求解分詞的狀態。 到這裏先引出一個比較簡單的問題，0號問題： 0.知道骰子的種類，骰子是什麼，每次扔什麼骰子，根據這個結果，求投出這個結果的機率是多少。 其實就是該知道的都知道了，直接求機率。

求解無非就是機率相乘了：

###③三個問題的簡單解答 ####1.看見了觀測序列，求隱藏序列 這裏解釋的就是第一種解法了，最大似然估計，也就是暴力解法。首先我知道我有3個骰子，六面骰子，八面骰子，四面骰子，同時結果我也知道了（1 6 3 5 2 7 3 5 2 4），如今想知道我每一次投的哪個骰子，是六面的仍是四面的仍是八面的呢？ 最簡單的方法就是窮舉了，從零個一直到第最後一個一次把每個機率算出來，鏈不長還行，鏈要是長了就算不了了，窮舉的數量太大。 接下來是要討論另外一個比較牛逼的算法，viterbi algorithm。 首先只扔一次骰子：

結果是1的話那麼機率最大的就是四面骰子了，1/4，其餘的都是1/6和1/8。接着再扔一次：

這個時候咱們就要計算三個值了，分別是D6，D8，D4的機率。要取到最大的機率，那麼第一次確定就是取到D4了，因此取到D6的最大機率：$$P_2(D6) = P(D4) P(D4->1)P(D6)P(D6->6) = \frac{1}{3}\frac{1}{4}\frac{1}{3}\frac{1}{6}$$ 取到D8的最大機率：$$P_2(D8) = P(D4) P(D4->1)P(D8)P(D8->6) = \frac{1}{3}\frac{1}{4}\frac{1}{3}\frac{1}{8}$$取到D4這輩子都不可能了，四面骰子沒有6。因此最大的序列就是D4，D6 繼續拓展一下，投多一個：

這個時候又要計算三種狀況了，分別計算爲D4，D6，D8的機率是多少：

能夠看到，機率最大的就是D4了，因此三次投骰子機率最大的隱含序列就是D4，D6，D4。因此不管這個序列多長，直接遞推計算能夠算出來的，算到最後一位的時候，直接反着找就行了。 ####2.求可觀測序列的機率 若是你的骰子有問題，要怎麼檢測出來？首先，能夠先算一下正常骰子投出一段序列的機率，再算算不正常的投出來序列的機率便可。若是小於正常的，那麼就有多是錯的了。

好比結果如上，咱們這個時候就不是要求隱含序列了，而是求出現這個觀測序列的總機率是多少。首先若是是隻投一個骰子：

這時候三種骰子都有可能：

再投一次：

一樣是要計算總的分數。

再投一次：

就是這樣按照套路計算一遍再比較結果便可。 ####3.只是知道觀測序列，想知道骰子是怎麼樣的？ 這個算法在後面會細講，Baum-welch算法。 ###④HMM的公式角度 下面正式開始講解，上面只是過一個印象。 ####1.馬爾科夫模型 要看隱馬爾科夫天然先動馬爾科夫是什麼。已知如今有N個有序的隨機變量，根據貝葉斯他們的聯合分佈能夠寫成條件連乘：

再繼續拆：$$p(x_1,x_2,x_3...x_N) = \prod_{i=1}^{n}p(x_i|x_{i-1},...,x_1)$$ 馬爾科夫鏈就是指，序列中的任何一個隨機變量在給定它的前一個變量的分佈於更早的變量無關。也就是說當前的變量只與前一個變量有關，與更早的變量是沒有關係的。用公式表示：

代進去：$$Q(\lambda,\hat\lambda) = \sum_Iln\pi_{i_1}P(O,I|\hat\lambda) + \sum_I(\sum_{t=1}^{T-1}lna_{i_ti_{t+1}})P(O,I|\hat\lambda)+\sum_I(\sum_{t=1}^{T}lnb_{i_to_t})P(O,I|\hat\lambda)$$ 爲何要分紅三個部分呢？由於求導一下就沒一邊了，開心。 求： 有一個條件：。既然是有條件了，拉格朗日乘子法就用上了。化簡一下上式：$$ \sum_Iln\pi_{i_1}P(O,I|\hat\lambda) = \sum_Iln\pi_{i_1}P(O,i_1=i|\hat\lambda) $$拉格朗日乘子法：

求導： 上面那幾條式子所有加起來： 因此 求A： 這裏的條件就是，仍是延用上面的方法，拉格朗日求導，其實就是改一下上面的公式就行了。因此有：

求B： 這裏的條件就是，因而：

這樣就是整個更新算法了。 ######有監督學習 若是已經有了一堆標註好的數據，那就沒有必要再去玩baum-welch算法了。這個算法就很簡單了，多的不說，直接上圖：

#####問題3：若是有了參數 和觀測序列，求 的機率最大的隱含序列 是什麼。

上面骰子的例子已經舉過了，再講一個下雨天的例子：

day	rain	sun
rain	0.7	0.3
sun	0.4	0.6

轉移機率如上

day	walk	shop	clean
rain	0.1	0.4	0.5
sun	0.6	0.3	0.1

發射機率入上 初始機率 已知模型參數和觀測三天行爲(walk,shop,clean)，求天氣。 ①首先是初始化：
②次日： ③最後一天： 接着就是回溯了，查看最後一天哪一個最大，倒着找，最後就是2,1,1了，也就是(sun,rain,rain)。 ###⑤Hidden Markov Model代碼實現 ####1.工具類的實現

class Tool(object):
    infinite = float(-2**31)
    @staticmethod
    def log_normalize(a):
        s = 0
        for x in a:
            s += x
        if s == 0:
            print('Normalize error,value equal zero')
            return
        s = np.log(s)
        for i in range(len(a)):
            if a[i] == 0:
                a[i] = Tool.infinite
            else:
                a[i] = np.log(a[i]) - s
        return a

    @staticmethod
    def log_sum(a):
        if not a:
            return Tool.infinite
        m = max(a)
        s = 0
        for t in a:
            s += np.exp(t - m)
        return m + np.log(s)


    @staticmethod
    def saveParameter(pi,A, B, catalog):
        np.savetxt(catalog + '/' + 'pi.txt', pi)
        np.savetxt(catalog + '/' + 'A.txt', A)
        np.savetxt(catalog + '/' + 'B.txt', B)
複製代碼

準備了幾個靜態方法，一個就是log標準化，也就是把全部的數組都歸一化操做，變成一個機率，log算加和，保存矩陣。全部的操做都是使用log表示，若是單單是隻用數值表示的話，由於涉及到不少的機率相乘，不少時候就變很小，根本檢測不出。而用log以後下限會大不少。 ####2.有監督算法的實現 其實就是根據上面的公式統計便可。

def supervised(filename):
    ''' The number of types of lastest variable is four,0B(begin)|1M(meddle)|2E(end)|3S(sigle) :param filename:learning fron this file :return: pi A B matrix '''
    pi = [0]*4
    a = [[0] * 4 for x in range(4)]
    b = [[0] * 65535 for x in range(4)]
    f = open(filename,encoding='UTF-8')
    data = f.read()[3:]
    f.close()
    tokens = data.split(' ')

    #start training
    last_q = 2
    old_process = 0
    allToken = len(tokens)
    print('schedule : ')
    for k, token in enumerate(tokens):
        process = float(k) /float(allToken)
        if process > old_process + 0.1:
            print('%.3f%%' % (process * 100))
            old_process = process
        token = token.strip()
        n = len(token)
        #empty we will choose another
        if n <= 0:
            continue
        #if just only one
        if n == 1:
            pi[3] += 1
            a[last_q][3] += 1
            b[3][ord(token[0])] += 1
            last_q = 3
            continue
        #if not
        pi[0] += 1
        pi[2] += 1
        pi[1] += (n-2)
        #transfer matrix
        a[last_q][0] += 1
        last_q = 2
        if n == 2:
            a[0][2] += 1
        else:
            a[0][1] += 1
            a[1][1] += (n-3)
            a[1][2] += 1
        #launch matrix
        b[0][ord(token[0])] += 1
        b[2][ord(token[n-1])] += 1
        for i in range(1, n-1):
            b[1][ord(token[i])] += 1
    pi = Tool.log_normalize(pi)
    for i in range(4):
        a[i] = Tool.log_normalize(a[i])
        b[i] = Tool.log_normalize(b[i])
    return pi, a, b
複製代碼

按照公式來，一個一個單詞判斷，若是是單個的那天然就是single，因此當的時候直接就是。初始狀態是由於一開始確定是結束以後才接着開始的，因此天然就是2，end。以後都是比較常規了。 ####3.viterbi算法的實現

def viterbi(pi, A, B, o):
    ''' viterbi algorithm :param pi:initial matrix :param A:transfer matrox :param B:launch matrix :param o:observation sequence :return:I '''
    T = len(o)
    delta = [[0 for i in range(4)] for t in range(T)]
    pre = [[0 for i in range(4)] for t in range(T)]
    for i in range(4):
        #first iteration
        delta[0][i] = pi[i] + B[i][ord(o[0])]
    for t in range(1, T):
        for i in range(4):
            delta[t][i] = delta[t-1][0] + A[0][i]
            for j in range(1, 4):
                vj = delta[t-1][j] + A[0][j]
                if delta[t][i] < vj:
                    delta[t][i] = vj
                    pre[t][i] = j
            delta[t][i] += B[i][ord(o[t])]
    decode = [-1 for t in range(T)]
    q = 0
    for i in range(1, 4):
        if delta[T-1][i] > delta[T-1][q]:
            q = i
    decode[T-1] = q
    for t in range(T-2, -1, -1):
        q = pre[t+1][q]
        decode[t] = q
    return decode
複製代碼

根據上面的例子來就好，先找到轉移機率最大的，迭代到最後找到機率最大的序列以後倒着回來找便可。最後就獲得一串編碼，而後使用這段編碼來進行劃分。

def segment(sentence, decode):
    N = len(sentence)
    i = 0
    while i < N:
        if decode[i] == 0 or decode[i] == 1:
            j = i+1
            while j < N:
                if decode[j] == 2:
                    break
                j += 1
            print(sentence[i:j+1],"|",end=' ')
            i = j+1
        elif decode[i] == 3 or decode[i] == 2:  # single
            print (sentence[i:i + 1], "|", end=' ')
            i += 1
        else:
            print ('Error:', i, decode[i] , end=' ')
            i += 1
複製代碼

這個就是根據編碼劃分句子的函數了。首先是經過有監督的學習獲得參數，而後用viterbi算法獲得編碼序列，再用segment函數作劃分便可。

if __name__ == '__main__':
    pi = np.loadtxt('supervisedParam/pi.txt')
    A = np.loadtxt('supervisedParam/A.txt')
    B = np.loadtxt('supervisedParam/B.txt')
    f = open("../Data/novel.txt" , encoding='UTF-8')
    data = f.read()[3:]
    f.close()
    decode = viterbi(pi, A, B, data)
    segment(data, decode)
複製代碼

執行代碼。

效果固然能夠了，畢竟是有監督。無監督的就沒有這麼秀了。 ####4.baum-welch算法的實現 參考上面三個公式，除了B有點難更新以外其餘的都很簡單。

def cal_alpha(pi, A, B, o, alpha):
    print('start calculating alpha......')
    for i in range(4):
        alpha[0][i] = pi[i] + B[i][ord(o[0])]
    T = len(o)
    temp = [0 for i in range(4)]
    del i
    for t in range(1, T):
        for i in range(4):
            for j in range(4):
                temp[j] = (alpha[t-1][j] + A[j][i])
            alpha[t][i] = Tool.log_sum(temp)
            alpha[t][i] += B[i][ord(o[t])]
    print('The calculation of alpha have been finished......')

def cal_beta(pi, A, B, o, beta):
    print('start calculating beta......')
    T = len(o)
    for i in range(4):
        beta[T-1][i] = 1
    temp = [0 for i in range(4)]
    del i
    for t in range(T-2, -1, -1):
        for i in range(4):
            beta[t][i] = 0
            for j in range(4):
                temp[j] = A[i][j] + B[j][ord(o[t + 1])] + beta[t + 1][j]
            beta[t][i] += Tool.log_sum(temp)
    print('The calculation of beta have been finished......')

def cal_gamma(alpha, beta, gamma):
    print('start calculating gamma......')
    for t in range(len(alpha)):
        for i in range(4):
            gamma[t][i] = alpha[t][i] + beta[t][i]
        s = Tool.log_sum(gamma[t])
        for i in range(4):
            gamma[t][i] -= s
    print('The calculation of gamma have been finished......')

def cal_kesi(alpha, beta, A, B, o, ksi):
    print('start calculating ksi......')
    T = len(o)
    temp = [0 for i in range(16)]
    for t in range(T - 1):
        k = 0
        for i in range(4):
            for j in range(4):
                ksi[t][i][j] = alpha[t][i] + A[i][j] + B[j][ord(o[t+1])] + beta[t+1][j]
                temp[k] = ksi[t][i][j]
                k += 1
        s = Tool.log_sum(temp)
        for i in range(4):
            for j in range(4):
                ksi[t][i][j] -= s
    print('The calculation of kesi have been finished......')

複製代碼

的更新按照公式便可。

def update(pi, A, B, alpha, beta, gamma, ksi, o):
    print('start updating......')
    T = len(o)
    for i in range(4):
        pi[i] = gamma[0][i]
    s1 = [0 for x in range(T-1)]
    s2 = [0 for x in range(T-1)]
    for i in range(4):
        for j in range(4):
            for t in range(T-1):
                s1[t] = ksi[t][i][j]
                s2[t] = gamma[t][i]
            A[i][j] = Tool.log_sum(s1) - Tool.log_sum(s2)
    s1 = [0 for x in range(T)]
    s2 = [0 for x in range(T)]
    for i in range(4):
        for k in range(65536):
            if k%5000 == 0:
                print(i, k)
            valid = 0
            for t in range(T):
                if ord(o[t]) == k:
                    s1[valid] = gamma[t][i]
                    valid += 1
                s2[t] = gamma[t][i]
            if valid == 0:
                B[i][k] = -Tool.log_sum(s2)
            else:
                B[i][k] = Tool.log_sum(s1[:valid]) - Tool.log_sum(s2)
    print('baum-welch algorithm have been finished......')
複製代碼

最後再封裝一把：

def baum_welch(pi, A, B, filename):
    f = open(filename , encoding='UTF-8')
    sentences = f.read()[3:]
    f.close()
    T = len(sentences)   # 觀測序列
    alpha = [[0 for i in range(4)] for t in range(T)]
    beta = [[0 for i in range(4)] for t in range(T)]
    gamma = [[0 for i in range(4)] for t in range(T)]
    ksi = [[[0 for j in range(4)] for i in range(4)] for t in range(T-1)]
    for time in range(1000):
        print('Time : ', time)
        sentence = sentences
        cal_alpha(pi, A, B, sentence, alpha)
        cal_beta(pi, A, B, sentence, beta)
        cal_gamma(alpha, beta, gamma)
        cal_kesi(alpha, beta, A, B, sentence, ksi)
        update(pi, A, B, alpha, beta, gamma, ksi, sentence)
        Tool.saveParameter(pi, A, B, 'unsupervisedParam')
        print('Save matrix successfully!')
複製代碼

初始化矩陣：

def inite():
    pi = [random.random() for x in range(4)]
    Tool.log_normalize(pi)
    A = [[random.random() for y in range(4)] for x in range(4)]
    A[0][0] = A[0][3] = A[1][0] = A[1][3] = A[2][1] = A[2][2] = A[3][1] = A[3][2] = 0
    B = [[random.random() for y in range(65536)] for x in range(4)]
    for i in range(4):
        A[i] = Tool.log_normalize(A[i])
        B[i] = Tool.log_normalize(B[i])
    return pi , A , B
複製代碼

最後跑一遍就OK了。 無監督算法使用的是EM框架，確定會存在局部最大的問題和初始值敏感，跑了56次，用的谷歌GPU跑，代碼沒有通過優化，跑的賊慢。

最後的效果也是慘不忍睹。 ###總結 機率圖模型大多都是圍繞着三個問題展開的，求觀測序列的機率最大，求隱含序列的機率最大，求參數。MEMM,RCF大多都會圍繞這幾個問題展開。求觀測序列的機率，暴力求解是爲了理解模型，前向後向算法纔是真正有用的；機率最大的隱含序列viterbi算法，動態規劃的思想最後回溯回來查找；求參數就是套用EM框架求解。 HMM是屬於生成模型，首先求出，轉移機率和表現機率直接建模，也就是對樣本密度建模，而後才進行推理預測。事實上有時候不少NLP問題是和先後相關，而不是隻是和前面的相關，HMM這裏明顯是隻和前面的隱變量有關，因此仍是存在侷限性的。對於優化和模型優缺點等寫了MEMM和RCF再一塊兒總結。

###最後附上GitHub代碼： github.com/GreenArrow2…