隐马尔可夫模型(HMM)+天气问题代码实现

 天气问题的代码实现：

state = ([0,1]) # states = ('Rainy' -> 0, 'Sunny' -> 1)
observations = np.array([0, 1, 2]) # observations = ('walk' -> 0 , 'shop' -> 1, 'clean' -> 2)
start_probability = np.array([0.6, 0.4])
transition_probability = np.array([[0.7, 0.3],
                                   [0.4, 0.6]])
emission_probability = np.array([[0.1, 0.4, 0.5],
                                 [0.6, 0.3, 0.1]])
P = np.zeros((3,2,2))
pp = np.zeros((3,2))
tmp = np.zeros(2)
day = np.zeros(3)

模拟隐马尔可夫思路：

主要是重现数理逻辑，除了需要注意kij区分外，没有其他问题。

for k in range (0, 3): # day
    for i in range (0, 2): # last weather
        for j in range (0, 2): # now weather
            if k == 0:
                P[k, i, j] = start_probability[ state[j] ] * emission_probability[ state[j], observations[k] ]
                pp[ k, state[j] ] = P[k, 1, state[j]]
            else:
                P[k, state[i], state[j]] = pp[ k-1, state[i] ] * transition_probability[ state[i], state[j] ] * emission_probability[ state[j], observations[k] ]
                pp[ k, state[j] ] = max( P[k, state[0], state[j]], P[k, state[1], state[j] ])

结果分析+输出：

写了两个函数来完成，递归回溯路径和输出。特别是Print()函数有关于比较大小的部分应有更好的写法，但我目前没有写出完善的其他方法。（Record_Weather()函数有个思路是利用python实现类似于C的结构体数组来进行）

def Record_Weather(k, j): # day, nowweather
    day[k] = j
    if k>0 and P[k, 0, j] > P[k, 1, j]:
        Record_Weather(k-1, 0)
    if k>0 and P[k, 0, j] < P[k, 1, j]:
        Record_Weather(k-1, 1)

def Print(day):
    if P[2, 0, 0] > P[2, 1, 0]:
        tmp[0] = P[2, 0, 0]
    else:
        tmp[0] = P[2, 1, 0]
    if P[2, 0, 1] > P[2, 1, 1]:
        tmp[1] = P[2, 0, 0]
    else:
        tmp[1] = P[2, 1, 0]
    if tmp[0]>tmp[1]:
        Record_Weather(2,0)
    else:
        Record_Weather(2,1)
    for i in range (3):
        if day[i] == 0:
            print("Rainy")
        else:
            print("Sunny")

写的过程中还是比较艰难的，因为开始有意识地训练自己去写一些可读性更强的代码，不能仅仅是重现数理逻辑。写代码期间也拓展了np.argsort()和用类实现结构体的知识，算是有提升。

感觉自己是有很多不好的代码习惯……希望能找到方法改掉。

全部代码：

import random
import numpy as np

state = ([0,1]) # states = ('Rainy' -> 0, 'Sunny' -> 1)
observations = np.array([0, 1, 2]) # observations = ('walk' -> 0 , 'shop' -> 1, 'clean' -> 2)
start_probability = np.array([0.6, 0.4])
transition_probability = np.array([[0.7, 0.3],
                                   [0.4, 0.6]])
emission_probability = np.array([[0.1, 0.4, 0.5],
                                 [0.6, 0.3, 0.1]])

def Record_Weather(k, j): # day, nowweather
    day[k] = j
    if k>0 and P[k, 0, j] > P[k, 1, j]:
        Record_Weather(k-1, 0)
    if k>0 and P[k, 0, j] < P[k, 1, j]:
        Record_Weather(k-1, 1)

def Print(day):
    if P[2, 0, 0] > P[2, 1, 0]:
        tmp[0] = P[2, 0, 0]
    else:
        tmp[0] = P[2, 1, 0]
    if P[2, 0, 1] > P[2, 1, 1]:
        tmp[1] = P[2, 0, 0]
    else:
        tmp[1] = P[2, 1, 0]
    if tmp[0]>tmp[1]:
        Record_Weather(2,0)
    else:
        Record_Weather(2,1)
    for i in range (3):
        if day[i] == 0:
            print("Rainy")
        else:
            print("Sunny")

P = np.zeros((3,2,2))
pp = np.zeros((3,2))
tmp = np.zeros(2)
day = np.zeros(3)

for k in range (0, 3): # day
    for i in range (0, 2): # last weather
        for j in range (0, 2): # now weather
            if k == 0:
                P[k, i, j] = start_probability[ state[j] ] * emission_probability[ state[j], observations[k] ]
                pp[ k, state[j] ] = P[k, 1, state[j]]
            else:
                P[k, state[i], state[j]] = pp[ k-1, state[i] ] * transition_probability[ state[i], state[j] ] * emission_probability[ state[j], observations[k] ]
                pp[ k, state[j] ] = max( P[k, state[0], state[j]], P[k, state[1], state[j] ])

Print(day)