用numpy实现隐马尔科夫模型的概率计算问题(前向算法)

实现用例为李航统计学习方法第177页例10.2

import numpy as np

'''
   用numpy实现隐马尔科夫模型的概率计算，即已知三元模型λ=（A,B,π）和观测序列O，
计算在模型λ下观测序列O出现的概率P(O|λ)
'''

'''一，定义模型的三要素'''

'''状态的集合'''
def status_gener(status_len:int):
    status_list = []
    for status_index in range(status_len):
        status_list.append(status_len)
    return status_list

'''观测的集合'''
obs_list = np.array(["红", "白", "红"])

'''时刻数'''
T = 3

'''初始取哪个盒子的概率分布'''
init_prob_dist = np.array([0.2,0.4,0.4]) 

'''状态转移概率分布'''
state_tran_prob = np.array([[0.5, 0.2, 0.3],
                                [0.3, 0.5, 0.2],
                                [0.2, 0.3, 0.5]])

'''观测概率分布'''
obs_prob_dist =  np.array([[0.5, 0.5],
                           [0.4, 0.6],
                           [0.7, 0.3]])

'''二，利用前向算法计算'''

'''初值计算'''
def init_value_calcu(obs_list, obs_prob_dist, init_prob_dist): # 3个形参分别为观测序列、观测概率分布、初始状态概率分布
    if obs_list[0] == "红":
        obs_prob_dist_col = obs_prob_dist[:,0]  # 取出每一时刻取出某个盒子后拿出红球的概率
    elif obs_list[0] == "白":
        obs_prob_dist_col = obs_prob_dist[:,1]  #  取出每一时刻取出某个盒子后拿出白球的概率
    init_value = init_prob_dist * obs_prob_dist_col # 初始取哪个盒子的概率分布*每一时刻取出某个盒子后拿出白或红球的概率     
    return init_value


'''递推计算'''
def recur_calcu(obs_list, state_tran_prob):
    # print(len(state_tran_prob))
    init_state_value = init_value_calcu(obs_list, obs_prob_dist, init_prob_dist)
    # print(init_state_value)
    a = 0
    time_num = len(state_tran_prob)
    for obs_id in obs_list[1:]:
        state = np.array([])
        for time_num_id in range(time_num):
            if a < 1:
                prob_value = np.cumsum(init_state_value * state_tran_prob[:, time_num_id], axis=0)[-1]
                state = np.append(state, prob_value)
            else:
                prob_value = np.cumsum(state_value * state_tran_prob[:, time_num_id], axis=0)[-1]
                state = np.append(state, prob_value)
        if obs_id == "红":
            state_value = state * obs_prob_dist[:, 0]
        elif obs_id == "白":
            state_value = state * obs_prob_dist[:, 1]
        a += 1
    final_value = np.cumsum(state_value, axis=0)[-1]
    return final_value
  
if __name__ == "__main__":
    print(recur_calcu(obs_list, state_tran_prob))

Done！