《统计学习方法》第十章，隐马尔科夫模型

▶ 隐马尔科夫模型的三个问题

● 代码

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

dataSize = 200
trainRatio = 0.3
epsilon = 1E-10
randomSeed = 109

def dataSplit(dataX, dataY, part):                                                  # 将数据集分割为训练集和测试集
    return dataX[:part], dataY[:part], dataX[part:], dataY[part:]

def normalCDF(x, μList, σList):
    return np.exp(-(x - μList)**2 / (2 * σList**2)) / (np.sqrt(2 * np.pi) * σList)

def targetIndex(x, xList):                      # 二分查找 xList 中大于 x 的最小索引    
    if x < xList[0]:
        return 0
    lp = 0
    rp = len(xList) - 1
    while lp < rp - 1:
        mp = (lp + rp) >> 1
        if(x <= xList[mp]):
            rp = mp
        else:
            lp = mp
    return rp

def createData(state, obs):                     # 创建数据，输入状态数、观测数，输出初始分布、状态转移分布、观测分布
    np.random.seed(randomSeed)
    π = np.random.rand(state)
    π /= np.sum(π)
    A = np.random.rand(state, state)    
    A = np.array( A / np.mat(np.sum(A,1)).T)    # 行归一化        
    B = np.random.rand(state, obs)    
    B = np.array(B / np.mat(np.sum(B,1)).T)    

    print("π = ", π)
    print("A = 
", A)
    print("B = 
", B)
    return π, A, B

def createString(π, A, B, length, count = 1):   # 创建状态列和观测列
    state, obs = np.shape(B)
    accπ = np.nancumsum(π)
    accA = np.cumsum(A, 1)
    accB = np.cumsum(B, 1)

    if count == 1:                                                                  # count 为 1 时生成的 string 只有一维，否则二维
        stateString = np.zeros(length, dtype = int)
        obsString = np.zeros(length, dtype = int)
        stateString[0] = targetIndex(np.random.rand(), accπ)                        # 随机选择初态
        for t in range(1, length):                                                  # 生成状态列
            stateString[t] = targetIndex(np.random.rand(), accA[stateString[t-1]])  
        for t in range(length):                                                     # 生成观测列
            obsString[t] = targetIndex(np.random.rand(), accB[stateString[t]])

    else:
        stateString = np.zeros([count, length], dtype = int)
        obsString = np.zeros([count, length], dtype = int)
        stateString[:,0] = [ targetIndex(np.random.rand(), accπ) for i in range(count) ]
        for t in range(1, length):
            stateString[:,t] = [ targetIndex(np.random.rand(), accA[stateString[i,t-1]]) for i in range(count) ]
        for t in range(length):
            obsString[:,t] = [ targetIndex(np.random.rand(), accB[stateString[i,t]]) for i in range(count) ]

    #print("state = 
", stateString)
    #print("obs = 
", obsString)
    return stateString, obsString

def pObsForwardSimplified(π, A, B, obsString):                  # 前向算法计算观测序列概率
    state, obs = np.shape(B)
    length = len(obsString)

    α = π * B[:,obsString[0]]
    for i in range(1, length):
        α = np.dot(α, A) * B[:,obsString[i]]                    # 等价代码 α = np.sum(α * A.T, 1) * B[:,obsString[i]]                            
    return np.sum(α)

def pObsForward(π, A, B, obsString):                            # 前向算法计算观测序列概率，并且输出前向概率矩阵
    state, obs = np.shape(B)
    length = len(obsString)

    α = np.zeros([length,state])
    α[0] = π * B[:,obsString[0]]
    for i in range(1, length):
        α[i] = np.dot(α[i-1], A) * B[:,obsString[i]]            # α[i] = np.sum(α[i-1], A.T, 1) * B[:,obsString[i]]
    return np.sum(α[-1]), α

def pObsBackwardSimplified(π, A, B, obsString):                 # 后向算法计算观测序列概率
    state, obs = np.shape(B)
    length = len(obsString)

    β = np.ones(state)
    for i in range(1, length)[::-1]:
        β = np.dot(A, β * B[:,obsString[i]])                    # β = np.sum(A * β * B[:,obsString[i]], 1)
    return np.sum(π * β * B[:,obsString[0]])

def pObsBackward(π, A, B, obsString):                           # 后向算法计算观测序列概率，并且输出后向概率矩阵
    state, obs = np.shape(B)
    length = len(obsString)

    β = np.zeros([length,state])
    β[-1] = np.ones(state)
    for i in range(1, length)[::-1]:
        β[i-1] = np.dot(A, β[i] * B[:,obsString[i]])            # β[i-1] = np.sum(A * β[i] * B[:,obsString[i]], 1)
    return np.sum(π * B[:,obsString[0]] * β[0]), β

def pObsMixedSimplified(A, B, obsString, α, β):                 # 用 α 和 β 来算观测序列概率
    return np.dot(α[0], np.dot(A, B[:,obsString[1]] * β[1]))    # np.sum(α[0] * np.sum(A * B[:,obsString[1]] * β[1], 1))

def pObsMixed(A, B, obsString, α, β):                           # 用 α 和 β 来算观测序列概率，这 length - 1 个结果都相同
    length = len(α)
    return np.array([ np.dot(α[t], np.dot(A, B[:,obsString[t+1]] * β[t+1])) for t in range(length-1) ])    
    #return np.array([ np.sum(α[t] * np.sum(A * B[:,obsString[t+1]] * β[t+1] , 1)) for t in range(length-1) ])

def pΓ(α, β):                                                   # 求各时刻隐变量处于各状态的概率
    t = α * β
    return np.array( t / np.mat(np.sum(t, 1)).T )               

def pΞ(A, B, obsString, α, β):                                  # t 时刻处于状态 i 而 t+1 时刻处于状态 j 的概率，t = 0 ~ length-2
    state, obs = np.shape(B)
    length = len(α)

    res = np.zeros([ length - 1, state, state ])
    for t in range(length - 1):
        res[t] = np.tile(α[t],[state,1]).T * A * B[:,obsString[t+1]] * β[t+1]
    return res / pObsMixedSimplified(A, B, obsString, α, β)

def supervisedLearn(dataState, dataObs, state, obs):            # 监督学习
    count, length = np.shape(dataState)
    π = np.zeros(state)
    A = np.zeros([state, state])
    B = np.zeros([state, obs])    
    
    for i in dataState[:,0]:
        π[i] +=1
    for i,j in zip(dataState[:,:-1].flatten(), dataState[:,1:].flatten()):
        A[i,j] +=1
    for i,j in zip(dataState.flatten(),dataObs.flatten()):
        B[i,j] +=1

    π /= count                                         
    A = np.array( A / np.mat(np.sum(A,1)).T)
    B = np.array( B / np.mat(np.sum(B,1)).T)
    return π, A, B

def baumWelchLearn(dataObs, state, obs):                        # Baum - WelchLearn 学习
    count, length = np.shape(dataObs)                           # 存在的问题：如果某个观测序列没有覆盖所有可能的观测结果，学习会产生 nan？
    π = np.random.rand(state)                                   # 选择初始值
    A = np.random.rand(state, state)
    B = np.random.rand(state, obs)
    π /= np.sum(π)
    A = np.array( A / np.mat(np.sum(A,1)).T)
    B = np.array( B / np.mat(np.sum(B,1)).T)

    for line in dataObs:                                        
        α = pObsForward(π, A, B, line)[1]
        β = pObsBackward(π, A, B, line)[1]
        ξ = pΞ(A, B, line, α, β)
        γ = pΓ(α, β)
        π = γ[0]
        A = np.array( np.sum(ξ, 0) / np.mat(np.sum(γ[:-1], 0)).T )
        for k in range(obs):
            B[:,k] = np.sum(γ[np.where(line==k)], 0)            
        B = np.array( B / np.mat(np.sum(γ, 0)).T )
    return π, A, B

def appropritePredict(x, π, A, B):                              # 近似预测，依次计算出α，β，γ，取 γ 每行最大值所在列号
    γ = pΓ(pObsForward(π, A, B, x)[1], pObsBackward(π, A, B, x)[1])
    return np.argmax(γ, 1)

def viterbiPredict(obsString, π, A, B):                         # Viterbi 预测，动态规划
    state, obs = np.shape(B)
    length = len(obsString)

    δ = π * B[:, obsString[0]]                                  # δ 记录以各状态为结尾的最大概率，ψ 记录该最大概率对应的上一节点号
    ψ = np.zeros([length, state])
    for t in range(1,length):
        temp = δ * A.T * B[:, obsString[t]]
        ψ[t-1] = np.argmax(temp,1)
        δ = np.max(temp,1)
        
    trace = np.zeros(length, dtype=np.int)
    trace[-1] = np.argmax(δ)
    for i in range(length - 1)[::-1]:
        trace[i] = ψ[i, trace[i+1]]
    
    return trace

def test(state, obs, length):                                   # 单次测试
    π, A, B = createData(state, obs)
    allState, allObs = createString(π, A, B, length, dataSize)           
    
    if allObs.ndim > 1:
        obsString = allObs[0]
    else:
        obsString = allObs
    '''                                                         
    π = np.array([0.2,0.4,0.4])                                 # 调试用数据
    A = np.array([[0.5,0.2,0.3],[0.3,0.5,0.2],[0.2,0.3,0.5]])
    B = np.array([[0.5,0.5],[0.4,0.6],[0.7,0.3]])    
    obsString = [0,1,0]
    '''    
    pOFS = pObsForwardSimplified(π, A, B, obsString)
    pOBS = pObsBackwardSimplified(π, A, B, obsString)
    pOF = pObsForward(π, A, B, obsString)
    pOB = pObsBackward(π, A, B, obsString)
    pOMS = pObsMixedSimplified(A, B, obsString, pOF[1], pOB[1])
    pOM  = pObsMixed(A, B, obsString, pOF[1], pOB[1])
    print("pOFS = %f, pOBS = %f
pOF  = %f, pOB  = %f
pOMS  = %f"%(pOFS, pOBS, pOF[0], pOB[0], pOMS))
    print("α = 
", pOF[1])
    print("β = 
", pOB[1])
    print("pOM = 
", pOM)
    γ = pΓ(pOF[1], pOB[1])
    print("γ = 
", γ)
    ξ = pΞ(A, B, obsString, pOF[1], pOB[1])
    print("ξ = 
", ξ)   

    print("
supervisedLearn:")
    para = supervisedLearn(allState, allObs, state, obs)    
    print("π = ", para[0])
    print("A = 
", para[1])
    print("B = 
", para[2])
    
    print("
baumWelchLearn:")
    para = baumWelchLearn(allObs, state, obs)    
    print("π = ", para[0])
    print("A = 
", para[1])
    print("B = 
", para[2])    

    myResult = [ appropritePredict(x, π, A, B) for x in allObs ]
    #print(myResult)
    errorRatio1 = np.sum( (np.array(myResult) != allState).astype(int) ) / (length * dataSize)
    
    myResult = [ viterbiPredict(x, π, A, B) for x in allObs ]
    #print(myResult)
    errorRatio2 = np.sum( (np.array(myResult) != allState).astype(int) ) / (length * dataSize)
    
    print("state = %d, obs = %d, length = %d
errorRatioAppropritePredict = %4f, errorRatioViterbiPredict = %4f"%(state, obs, length, errorRatio1, errorRatio2))

if __name__ == '__main__':
    test(4, 3, 10)

● 输出结果1，成功复现树上的样例数据

pOFS = 0.130218, pOBS = 0.130218
pOF  = 0.130218, pOB  = 0.130218
pOMS  = 0.130218
α =
 [[0.1      0.16     0.28    ]
 [0.077    0.1104   0.0606  ]
 [0.04187  0.035512 0.052836]]
β =
 [[0.2451 0.2622 0.2277]
 [0.54   0.49   0.57  ]
 [1.     1.     1.    ]]
pOM =
 [0.130218 0.130218]
γ =
 [[0.18822283 0.32216744 0.48960973]
 [0.31931069 0.41542644 0.26526287]
 [0.32153773 0.27271191 0.40575036]]
ξ =
 [[[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]]]
appropritePredict: [2 1 2]
ViterbiPredict: [2 2 2]

● 输出结果2，用自己的数据来跑有问题【坑】

π =  [0.09536369 0.4423878  0.33021783 0.13203068]
A =
 [[0.38752484 0.17316366 0.39198697 0.04732453]
 [0.09835982 0.23793441 0.38519171 0.27851407]
 [0.44365676 0.12177898 0.3055296  0.12903465]
 [0.29614541 0.30691252 0.04969234 0.34724973]]
B =
 [[0.12418479 0.76243987 0.11337534]
 [0.34515219 0.28857853 0.36626928]
 [0.27430006 0.14931117 0.57638877]
 [0.22205449 0.24865406 0.52929144]]

supervisedLearn:
π =  [0.1  0.42 0.33 0.15]
A =
 [[0.38233515 0.18269812 0.38747731 0.04748941]
 [0.10843373 0.23782085 0.37139864 0.28234678]
 [0.44822934 0.12310287 0.30084317 0.12782462]
 [0.30813953 0.27151163 0.04476744 0.3755814 ]]
B =
 [[0.13203593 0.75479042 0.11317365]
 [0.34575569 0.313147   0.34109731]
 [0.27606952 0.14639037 0.57754011]
 [0.23502304 0.24193548 0.52304147]]

baumWelchLearn:
π =  [1.60211593e-01 8.26662258e-01 8.91008964e-06 1.31172391e-02]
A =
 [[0.17940016 0.25527889 0.31784884 0.24747211]
 [0.39215007 0.15576319 0.22874897 0.22333777]
 [0.32884167 0.30016482 0.15265056 0.21834295]
 [0.36968335 0.03746577 0.34497214 0.24787874]]
B =
 [[0.34244911 0.24148873 0.41606216]
 [0.29723256 0.24124855 0.46151888]
 [0.18049919 0.57465181 0.244849  ]
 [0.29483333 0.37769373 0.32747294]]
state = 4, obs = 3, length = 100
errorRatioAppropritePredict = 0.489400, errorRatioViterbiPredict = 0.692400