隐马尔科夫模型

# coding=utf8import sysimport osimport numpy as np#region 运行环境reload(sys)sys.setdefaultencoding('utf-8')#endregiondef forward(A, B, PI, O):    """    算法10.2 观测序列概率前向算法    :param A: 状态转移概率矩阵 N*N    :param B: 观测概率矩阵 N*M    :param PI: 初始化状态概率 N    :param O: 观测序列 T    :return: P(0|A,B,pi)    """    N=A.shape[0]#可能的状态值数目    T = len(O)  # 时间序列长度    alpha_list=np.zeros([N,T])    alpha_list[:,0]= PI * B[:, O[0] - 1]    for t in range(2,T+1):        alpha_list[:,t-1]=np.matmul(alpha_list[:,t-2],A)*B[:,O[t-1]-1]    p=np.sum(alpha_list[:,T-1])    return  alpha_list,pdef backard(A,B,PI,O):    """    算法10.3 贯彻序列概率的后向算法    :param A: 状态转移概率矩阵 N*N    :param B: 观测概率矩阵 N*M    :param PI: 初始化状态概率 N    :param O: 观测序列 T    :return: P(0|A,B,pi)    """    N=A.shape[0]#可能的状态值数目    T=len(O)#时间序列长度    beta_list=np.zeros([N,T])    beta_list[:,T-1]=np.ones([1,N])    for t in range(T,1,-1):        beta_list[:,t-2]=np.matmul(B[:,O[t-1]-1]*beta_list[:,t-1],A.T)    p=np.sum(PI*B[:,O[0]-1]*beta_list[:,0])    return beta_list,pdef initP(M, N):    """    初始化一个M*N的概率矩阵    :param M:    :param N:    :return:    """    if M==0 or N==0:        return None    max_dim=M if M>N else N    A=np.random.uniform(0,1.0/max_dim,[M,N])    # A[M - 1, :] = 1 - np.sum(A[0:M - 1, :], axis=0)    # A[:, N - 1] = 1 - np.sum(A[:, 0:N - 1], axis=1)    return A# def calQ(PI,A,B,I,O):#     n,t=I.shape#     Q=0#     for i in range(n):#         p=PI[I[i,0]-1]*B[I[i,0]-1,O[i,0]-1]#         asum=0#         bsum=0#         for j in range(1,t):#             p=p*A[I[i,t-1]-1,I[i,t]-1]*B[I[i,0]-1,O[i,0]-1]#             asum+=np.log(A[I[i,t-1]-1,I[i,t]-1])#             bsum+=np.log(B[I[i,0]-1,O[i,0]-1])#         Q=Q+(np.log(PI[I[i,0]-1])+asum+bsum)*p#     return Qdef calGamma(A, B, PI, O):    """    计算概率P(O,i(t)=q(i)|A,B,PI)的N*T矩阵    :param A:    :param B:    :param PI:    :param O:    :return:    """    alpha,p = forward(A, B, PI, O)    beta,_= backard(A, B, PI, O)    N=A.shape[0]    T=len(O)    gamma=np.zeros([N,T])    for i in range(N):        for t in range(T):            gamma[i,t]=alpha[i,t]*beta[i,t]/p    return gammadef calEPS(A, B, PI, O):    """    计算概率P(O,i(t)=q(i),i(t+1)=q(i+1)|A,B,PI)的N*N*T矩阵    :param A:    :param B:    :param PI:    :param O:    :return:    """    alpha,p = forward(A, B, PI, O)    beta,_= backard(A, B, PI, O)    N,T=alpha.shape    eps=np.zeros([T-1,N,N])    for t in range(T - 1):        for i in range(N):            for j in range(N):                eps[t,i,j]=alpha[i,t]*A[i,j]*B[j,O[t+1]-1]*beta[j,t+1]/p    return epsdef baum_welch(O,M,N,precision=1e-6):    """    算法10.4 隐马尔科夫模型参数学习    :param O:观测序列    :param M:观测值数    :param N:状态数    :param precision:精确度    :return:A,B,PI    """    T=len(O)    A=initP(N, N)    B = initP(N, M)    PI = initP(1, N).flatten()    A_new=A.copy()    B_new=B.copy()    PI_new=PI.copy()    n=0    while True:        gamma=calGamma(A, B, PI, O)        eps=calEPS(A, B, PI, O)        for i in range(N):            PI_new[i]=gamma[i,0]        for i in range(N):            for j in range(N):                A_new[i,j]=np.sum(eps[:,i,j])/np.sum(gamma[i,0:T-1])        for i in range(N):            for k in range(M):                B_new[i,k]=np.sum(gamma[i,:][O==k+1])/np.sum(gamma[i,:])        #收敛判断        PI_add=np.linalg.norm(PI_new-PI)        A_add=np.linalg.norm(A_new-A)        B_add = np.linalg.norm(B_new-B)        if PI_add<precision and A_add<precision and B_add<precision:            break        A=A_new.copy()#注意要使用拷贝        B=B_new.copy()        PI=PI_new.copy()        n=n+1    print(n)    return A,B,PIdef viterbi(A,B,PI,O):    """    算法10.5 维特比算法做预测    :param A:状态转移概率矩阵    :param B:观测概率矩阵    :param PI:初始化状态概率矩阵    :param O:观测序列    :return:I 状态序列    """    N=A.shape[0]    M=B.shape[1]    T=len(O)    deltaArr=np.zeros([T,N])    psiArr=np.zeros([T,N],dtype=np.int)    deltaArr[0,:]=PI*B[:,O[0]-1]    for t in range(1,T):        for i in range(N):            deltaArr[t,i]=np.max(deltaArr[t-1,:]*A[:,i])*B[i,O[t]-1]            psiArr[t,i]=np.argmax(deltaArr[t-1,:]*A[:,i])    p=np.max(deltaArr[T-1,:])    it=np.argmax(deltaArr[T-1,:])    I=np.zeros(T,dtype=np.int)    I[T-1]=it    for t in range(T-1,0,-1):        it=psiArr[t,it]        I[t-1]=it    I=I+1    return Idef main():    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]])    PI = np.array([0.2, 0.4, 0.4])    O = np.array([1, 2,1])    # alpha = forward(A, B, PI, O)    # print(alpha)    # beta = backard(A, B, PI, O)    # print(beta)    # gamma=calGamma(A, B, PI, O)    # print(gamma)    # eps=calEPS(A, B, PI, O)    # print(eps)    # N=3    # M=2    # A_new,B_new,PI_new=baum_welch(O,M,N)    # print(A_new)    # print(B_new)    # print(PI_new)    I=viterbi(A,B,PI,O)    print(I)main()