HHM(Hidden Markov Model)-隐式马尔科夫模型

HHM(Hidden Markov Model)-隐式马尔科夫模型

隐式马尔科夫模型是通过对观察序列的观测推导出隐藏状态。举例说明，由于本人是在十一小长假写下该篇文章，因此就拿小长假的出行计划介绍该算法，通常本人会根据特定的天气采取不同的出行计划，如下图所示。在这7天的小长假中，2天阴天（1天打扫房间，1天逛商场），3天晴天（2天打篮球，1天逛商场），2天雨天（2天全部打扫房间），假设现在我们无法获取天气情况（可以理解成你被关进了小黑屋），但是你可以获取到我的行为数据，现在让你根据行为数据推导出天气情况，这个过程就是马尔可夫模型，其中行为数据就是观察序列，天气就是隐藏状态，这个推导过程是根据历史的数据，通过统计的方式计算出各种概率，通过概率的方式根据观察序列推导隐藏状态。

整个隐式马尔科夫模型的框架如下图所示，通过对观察序列（时间序列）的观察，推导出隐藏状态。这里需要注意一点，观察序列一定是受到隐藏状态的影响（箭头从隐藏状态指向了观察序列），也就是要有一定的因（隐藏状态）果（观察序列）关系，需要注意的是这种因果关系要有一定的时效性，不能说今天的天气影响到我后天的行为，所以当使用马尔可夫模型时，一定要注意这一点。

算法讲解:

上文提到HMM是通过概率统计来进行推导的，整个推导过程可以用五个元素来描述，即2个状态集合和3个概率矩阵：

1. 隐藏状态S：

之所以叫隐藏状态，是因为无法直接观测到，只能通过历史发生后才能知道，就好比上面的天气情况，这里的天气是无法直接观测到的（假设关进小黑屋）。

S{阴天，晴天，雨天}

2. 可观测状态O：

在模型中与隐藏状态相关联，与隐藏状态的关系是因果关系，可以通过观测而得到。

O{打扫房间，打篮球，逛商场}

3. 隐藏状态初始概率矩阵兀：

通过对历史数据的分析，可以统计出每个隐藏状态的概率，就比如上面的十一小长假（2天阴天/7，3天晴天/7，2天雨天/7）：

4. 隐藏状态转移概率矩阵A:

这里统计的是隐藏状态之间的转移概率，即一个状态转移到其他状态的概率，假设下一个状态只和上一个状态相关，就是一阶马尔科夫链，和前n个状态相关，就是n阶马尔科夫链，HMM通常假设是一阶马尔可夫链。由于我们只有三个隐藏状态，因此是一个3X3的二维矩阵：

十一假期天气分布：

转移概率统计：

5. 发射概率B：

这里统计的是观测值在所有隐藏状态下的概率，是将观察序列与隐藏状态进行关联统计，相当于统计因果关系概率，即在某个状态下，所表现出来的观测值的概率，只与当前时刻的状态有关，与其他时候的状态和观测值无关。

发射概率的公式以及统计如下，在某个隐藏状态Sj下，所有观测序列的概率之和等于1，qt来自隐藏状态的值，vk来自观察序列。

当然，也有的开发团队会利用贝叶斯理论，乘上先验概率：

根据上面统计的结果，可以利用viterbi算法求出针对观测序列（obs）所对应的最可能的隐藏状态（states），框架如下图所示。

需要注意的是，需要将上面的概率进行ln处理，这样可以将概率相乘转成概率相加。

viterbi算法简单来说就是求出最短路径，这个最短路径代表了最大概率(概率公式如下)，通过最短路径（最大概率）求出观察序列所对应的隐藏状态，该算法是将初始概率，转移概率和发射概率相结合算出的最大概率，具体算法如下代码所示。

概率公式：

框架：

viterbi代码示例：

obs:观测序列

states：隐藏状态

start_p: 隐藏状态初始概率

trans_p:隐藏状态转移概率

emit_p: 发射概率

PrevStatus：上一时刻可能的隐藏状态

def viterbi(obs, states, start_p, trans_p, emit_p):
  V = [{}]  
  path = {}
  for y in states:  
    V[0][y] = start_p[y] + emit_p[y].get(obs[0], MIN_FLOAT)
    path[y] = [y]
  for t in xrange(1, len(obs)):
    V.append({})
    newpath = {}
    for y in states:
      em_p = emit_p[y].get(obs[t], MIN_FLOAT)
      (prob, state) = max(
        [(V[t - 1][y0] + trans_p[y0].get(y, MIN_FLOAT) + em_p, y0) for y0 in PrevStatus[y]])
      V[t][y] = prob
      newpath[y] = path[state] + [y]
    path = newpath
  (prob, state) = max((V[len(obs) - 1][y], y) for y in 'ES')
  return (prob, path[state])

测试代码：

from hmmlearn import hmm
import numpy as np
states = ["Rainy", "Sunny", "Cloudy"]
n_states = len(states)
observations = ["basketball", "shop", "clean"]
n_observations = len(observations)
start_probability = np.array([0.28, 0.44, 0.28])
transition_probability = np.array([[0, 1, 0 ],
                                  [0, 0.5,0.5],
                                  [1, 0, 0]])
emission_probability = np.array([
  [0, 0, 1],
  [0.67, 0.33, 0],
  [0, 0.5, 0.5]
  ])        
model = hmm.MultinomialHMM(n_components=n_states)
model.n_features = len(observations)
model.startprob_= start_probability
model.transmat_ = transition_probability
model.emissionprob_ = emission_probability

#给出一个可见序列,0=basketball, 1=shop, 2= clean#logprob#该参数反映模型拟合的好坏,数值越大越好jack_Actions = np.array([[2, 0, 1, 1, 2, 0]]).T
logprob, weathers = model.decode(jack_Actions , algorithm="viterbi")
print("Jack Actions:", ", ".join(map(lambda x: observations[x[0]], jack_Actions)))
print("weathers:", ", ".join(map(lambda x: states[x], weathers)))