python模擬隱馬爾可夫模型的方法是什么

這篇文章主要講解了“python模擬隱馬爾可夫模型的方法是什么”，文中的講解內(nèi)容簡(jiǎn)單清晰，易于學(xué)習(xí)與理解，下面請(qǐng)大家跟著小編的思路慢慢深入，一起來研究和學(xué)習(xí)“python模擬隱馬爾可夫模型的方法是什么”吧！

import numpy as np


class HiddenMarkov:
    def __init__(self):
        self.alphas = None
        self.forward_P = None
        self.betas = None
        self.backward_P = None

    # 前向算法
    # Q 是狀態(tài)集合，里面包含了所有可能的狀態(tài)
    # V 是我們的觀測(cè)的集合，里面包含了所有可能的觀測(cè)結(jié)果
    # A 狀態(tài)轉(zhuǎn)移概率分布
    # B 觀測(cè)概率分布
    # O 觀測(cè)序列，依次為觀測(cè)值
    # PI 初始概率分布。根據(jù)這個(gè)先生成初始狀態(tài)。
    def forward(self, Q, V, A, B, O, PI):
        # 狀態(tài)序列的大小
        N = len(Q)
        # 觀測(cè)序列的大小
        M = len(O)
        # 初始化前向概率alpha值
        alphas = np.zeros((N, M))
        # 時(shí)刻數(shù)=觀測(cè)序列數(shù)
        T = M
        # 遍歷每一個(gè)時(shí)刻，計(jì)算前向概率alpha值
        for t in range(T):
            # 得到序列對(duì)應(yīng)的索引
            indexOfO = V.index(O[t])
            # 遍歷狀態(tài)序列
            for i in range(N):
                # 初始化alpha初值
                if t == 0:
                    # P176 公式(10.15)
                    alphas[i][t] = PI[t][i] * B[i][indexOfO]
                    print('alpha1(%d) = p%db%db(o1) = %f' %
                          (i + 1, i, i, alphas[i][t]))
                else:
                    # P176 公式(10.16)
                    alphas[i][t] = np.dot([alpha[t - 1] for alpha in alphas],
                                          [a[i] for a in A]) * B[i][indexOfO]
                    print('alpha%d(%d) = [sigma alpha%d(i)ai%d]b%d(o%d) = %f' %
                          (t + 1, i + 1, t - 1, i, i, t, alphas[i][t]))
        # P176 公式(10.17)
        self.forward_P = np.sum([alpha[M - 1] for alpha in alphas])
        self.alphas = alphas

    # 后向算法
    # Q 是狀態(tài)集合，里面包含了所有可能的狀態(tài)
    # V 是我們的觀測(cè)的集合，里面包含了所有可能的觀測(cè)結(jié)果
    # A 狀態(tài)轉(zhuǎn)移概率分布
    # B 觀測(cè)概率分布
    # O 觀測(cè)序列，依次為觀測(cè)值
    # PI 初始概率分布。根據(jù)這個(gè)先生成初始狀態(tài)。
    def backward(self, Q, V, A, B, O, PI):
        # 狀態(tài)序列的大小
        N = len(Q)
        # 觀測(cè)序列的大小
        M = len(O)
        # 初始化后向概率beta值，P178 公式(10.19)
        betas = np.ones((N, M))
        for i in range(N):
            print('beta%d(%d) = 1' % (M, i + 1))
        # 對(duì)觀測(cè)序列逆向遍歷
        for t in range(M - 2, -1, -1):
            # 得到序列對(duì)應(yīng)的索引
            indexOfO = V.index(O[t + 1])
            # 遍歷狀態(tài)序列
            for i in range(N):
                # P178 公式(10.20)
                betas[i][t] = np.dot(
                    np.multiply(A[i], [b[indexOfO] for b in B]),
                    [beta[t + 1] for beta in betas])
                realT = t + 1
                realI = i + 1
                print('beta%d(%d) = sigma[a%djbj(o%d)beta%d(j)] = (' %
                      (realT, realI, realI, realT + 1, realT + 1),
                      end='')
                for j in range(N):
                    print("%.2f * %.2f * %.2f + " %
                          (A[i][j], B[j][indexOfO], betas[j][t + 1]),
                          end='')
                print("0) = %.3f" % betas[i][t])
        # 取出第一個(gè)值
        indexOfO = V.index(O[0])
        self.betas = betas
        # P178 公式(10.21)
        P = np.dot(np.multiply(PI, [b[indexOfO] for b in B]),
                   [beta[0] for beta in betas])
        self.backward_P = P
        print("P(O|lambda) = ", end="")
        for i in range(N):
            print("%.1f * %.1f * %.5f + " %
                  (PI[0][i], B[i][indexOfO], betas[i][0]),
                  end="")
        print("0 = %f" % P)

    # 維特比算法：動(dòng)態(tài)規(guī)劃解隱馬爾代夫模型預(yù)測(cè)問題
    # Q 是狀態(tài)集合，里面包含了所有可能的狀態(tài)
    # V 是我們的觀測(cè)的集合，里面包含了所有可能的觀測(cè)結(jié)果
    # A 狀態(tài)轉(zhuǎn)移概率分布
    # B 觀測(cè)概率分布
    # O 觀測(cè)序列，依次為觀測(cè)值
    # PI 初始概率分布。根據(jù)這個(gè)先生成初始狀態(tài)。
    def viterbi(self, Q, V, A, B, O, PI):
        # 狀態(tài)序列的大小
        N = len(Q)
        # 觀測(cè)序列的大小
        M = len(O)
        # 初始化daltas：存當(dāng)前時(shí)刻當(dāng)前狀態(tài)的所有單個(gè)路徑的概率最大值
        deltas = np.zeros((N, M))
        # 初始化psis：存當(dāng)前時(shí)刻當(dāng)前狀態(tài)所有單個(gè)路徑中概率最大路徑的前一時(shí)刻結(jié)點(diǎn)
        psis = np.zeros((N, M))
        # 初始化最優(yōu)路徑矩陣，該矩陣維度與觀測(cè)序列維度相同。這是我們最后的輸出。
        I = np.zeros((1, M))
        # 遍歷觀測(cè)序列
        for t in range(M):
            # 遞推從t=2開始
            realT = t + 1
            # 得到序列對(duì)應(yīng)的索引
            indexOfO = V.index(O[t])
            for i in range(N):
                realI = i + 1
                if t == 0:
                    # P185 算法10.5 步驟(1)
                    deltas[i][t] = PI[0][i] * B[i][indexOfO]
                    psis[i][t] = 0
                    print('delta1(%d) = pi%d * b%d(o1) = %.2f * %.2f = %.2f' %
                          (realI, realI, realI, PI[0][i], B[i][indexOfO],
                           deltas[i][t]))
                    print('psis1(%d) = 0' % (realI))
                else:
                    # # P185 算法10.5 步驟(2)
                    deltas[i][t] = np.max(
                        np.multiply([delta[t - 1] for delta in deltas],
                                    [a[i] for a in A])) * B[i][indexOfO]
                    print(
                        'delta%d(%d) = max[delta%d(j)aj%d]b%d(o%d) = %.2f * %.2f = %.5f'
                        % (realT, realI, realT - 1, realI, realI, realT,
                           np.max(
                               np.multiply([delta[t - 1] for delta in deltas],
                                           [a[i] for a in A])), B[i][indexOfO],
                           deltas[i][t]))
                    # 對(duì)于y=f(x)，argmax返回取得最大值y時(shí)的x
                    psis[i][t] = np.argmax(
                        np.multiply([delta[t - 1] for delta in deltas],
                                    [a[i] for a in A]))
                    print('psis%d(%d) = argmax[delta%d(j)aj%d] = %d' %
                          (realT, realI, realT - 1, realI, psis[i][t]))
        # 得到最優(yōu)路徑的終結(jié)點(diǎn)
        I[0][M - 1] = np.argmax([delta[M - 1] for delta in deltas])
        print('i%d = argmax[deltaT(i)] = %d' % (M, I[0][M - 1] + 1))
        # 遞歸由后向前得到其他結(jié)點(diǎn)
        for t in range(M - 2, -1, -1):
            I[0][t] = psis[int(I[0][t + 1])][t + 1]
            print('i%d = psis%d(i%d) = %d' %
                  (t + 1, t + 2, t + 2, I[0][t] + 1))
        # 輸出最優(yōu)路徑
        print('最優(yōu)路徑是：', "->".join([str(int(i + 1)) for i in I[0]]))


# 習(xí)題10.1
Q = [1, 2, 3]
V = ['紅', '白']
A = [[0.5, 0.2, 0.3], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]]
B = [[0.5, 0.5], [0.4, 0.6], [0.7, 0.3]]
# O = ['紅', '白', '紅', '紅', '白', '紅', '白', '白']
O = ['紅', '白', '紅', '白']    # 習(xí)題10.1的例子
PI = [[0.2, 0.4, 0.4]]
HMM = HiddenMarkov()
HMM.forward(Q, V, A, B, O, PI)
print("P(O|λ)={}".format(HMM.forward_P))
# HMM.backward(Q, V, A, B, O, PI)
# print("P(O|λ)={}".format(HMM.backward_P))
# HMM.viterbi(Q, V, A, B, O, PI)

• 結(jié)果

1. 前向算法

alpha1(1) = p0b0b(o1) = 0.100000
alpha1(2) = p1b1b(o1) = 0.160000
alpha1(3) = p2b2b(o1) = 0.280000
alpha2(1) = [sigma alpha0(i)ai0]b0(o1) = 0.077000
alpha2(2) = [sigma alpha0(i)ai1]b1(o1) = 0.110400
alpha2(3) = [sigma alpha0(i)ai2]b2(o1) = 0.060600
alpha3(1) = [sigma alpha1(i)ai0]b0(o2) = 0.041870
alpha3(2) = [sigma alpha1(i)ai1]b1(o2) = 0.035512
alpha3(3) = [sigma alpha1(i)ai2]b2(o2) = 0.052836
alpha4(1) = [sigma alpha2(i)ai0]b0(o3) = 0.021078
alpha4(2) = [sigma alpha2(i)ai1]b1(o3) = 0.025188
alpha4(3) = [sigma alpha2(i)ai2]b2(o3) = 0.013824
P(O|λ)=0.06009079999999999

2. 后向算法

beta4(1) = 1
beta4(2) = 1
beta4(3) = 1
beta3(1) = sigma[a1jbj(o4)beta4(j)] = (0.50 * 0.50 * 1.00 + 0.20 * 0.60 * 1.00 + 0.30 * 0.30 * 1.00 + 0) = 0.460
beta3(2) = sigma[a2jbj(o4)beta4(j)] = (0.30 * 0.50 * 1.00 + 0.50 * 0.60 * 1.00 + 0.20 * 0.30 * 1.00 + 0) = 0.510
beta3(3) = sigma[a3jbj(o4)beta4(j)] = (0.20 * 0.50 * 1.00 + 0.30 * 0.60 * 1.00 + 0.50 * 0.30 * 1.00 + 0) = 0.430
beta2(1) = sigma[a1jbj(o3)beta3(j)] = (0.50 * 0.50 * 0.46 + 0.20 * 0.40 * 0.51 + 0.30 * 0.70 * 0.43 + 0) = 0.246
beta2(2) = sigma[a2jbj(o3)beta3(j)] = (0.30 * 0.50 * 0.46 + 0.50 * 0.40 * 0.51 + 0.20 * 0.70 * 0.43 + 0) = 0.231
beta2(3) = sigma[a3jbj(o3)beta3(j)] = (0.20 * 0.50 * 0.46 + 0.30 * 0.40 * 0.51 + 0.50 * 0.70 * 0.43 + 0) = 0.258
beta1(1) = sigma[a1jbj(o2)beta2(j)] = (0.50 * 0.50 * 0.25 + 0.20 * 0.60 * 0.23 + 0.30 * 0.30 * 0.26 + 0) = 0.112
beta1(2) = sigma[a2jbj(o2)beta2(j)] = (0.30 * 0.50 * 0.25 + 0.50 * 0.60 * 0.23 + 0.20 * 0.30 * 0.26 + 0) = 0.122
beta1(3) = sigma[a3jbj(o2)beta2(j)] = (0.20 * 0.50 * 0.25 + 0.30 * 0.60 * 0.23 + 0.50 * 0.30 * 0.26 + 0) = 0.105
P(O|lambda) = 0.2 * 0.5 * 0.11246 + 0.4 * 0.4 * 0.12174 + 0.4 * 0.7 * 0.10488 + 0 = 0.060091
P(O|λ)=[0.0600908]

3. 維特比算法

delta1(1) = pi1 * b1(o1) = 0.20 * 0.50 = 0.10
psis1(1) = 0
delta1(2) = pi2 * b2(o1) = 0.40 * 0.40 = 0.16
psis1(2) = 0
delta1(3) = pi3 * b3(o1) = 0.40 * 0.70 = 0.28
psis1(3) = 0
delta2(1) = max[delta1(j)aj1]b1(o2) = 0.06 * 0.50 = 0.02800
psis2(1) = argmax[delta1(j)aj1] = 2
delta2(2) = max[delta1(j)aj2]b2(o2) = 0.08 * 0.60 = 0.05040
psis2(2) = argmax[delta1(j)aj2] = 2
delta2(3) = max[delta1(j)aj3]b3(o2) = 0.14 * 0.30 = 0.04200
psis2(3) = argmax[delta1(j)aj3] = 2
delta3(1) = max[delta2(j)aj1]b1(o3) = 0.02 * 0.50 = 0.00756
psis3(1) = argmax[delta2(j)aj1] = 1
delta3(2) = max[delta2(j)aj2]b2(o3) = 0.03 * 0.40 = 0.01008
psis3(2) = argmax[delta2(j)aj2] = 1
delta3(3) = max[delta2(j)aj3]b3(o3) = 0.02 * 0.70 = 0.01470
psis3(3) = argmax[delta2(j)aj3] = 2
delta4(1) = max[delta3(j)aj1]b1(o4) = 0.00 * 0.50 = 0.00189
psis4(1) = argmax[delta3(j)aj1] = 0
delta4(2) = max[delta3(j)aj2]b2(o4) = 0.01 * 0.60 = 0.00302
psis4(2) = argmax[delta3(j)aj2] = 1
delta4(3) = max[delta3(j)aj3]b3(o4) = 0.01 * 0.30 = 0.00220
psis4(3) = argmax[delta3(j)aj3] = 2
i4 = argmax[deltaT(i)] = 2
i3 = psis4(i4) = 2
i2 = psis3(i3) = 2
i1 = psis2(i2) = 3
最優(yōu)路徑是： 3->2->2->2

感謝各位的閱讀，以上就是“python模擬隱馬爾可夫模型的方法是什么”的內(nèi)容了，經(jīng)過本文的學(xué)習(xí)后，相信大家對(duì)python模擬隱馬爾可夫模型的方法是什么這一問題有了更深刻的體會(huì)，具體使用情況還需要大家實(shí)踐驗(yàn)證。這里是億速云，小編將為大家推送更多相關(guān)知識(shí)點(diǎn)的文章，歡迎關(guān)注！