隐马尔科夫模型

最新推荐文章于 2024-08-12 22:49:39 发布

HumorousJack

最新推荐文章于 2024-08-12 22:49:39 发布

阅读量150

点赞数

分类专栏：机器学习文章标签： python 机器学习源码隐马尔科夫模型

本文链接：https://blog.csdn.net/kylinat2688/article/details/78668488

版权

机器学习专栏收录该内容

5 篇文章 0 订阅

订阅专栏

# coding=utf8
import sys
import os
import numpy as np

#region 运行环境
reload(sys)
sys.setdefaultencoding('utf-8')
#endregion

def 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,p

def 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,p

def 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 Q

def 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 gamma

def 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 eps

def 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,PI

def 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 I

def 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()