HMM

https://zhuanlan.zhihu.com/p/62529131

一、知识回顾：

HMM(隐马尔可夫模型)的完整参数集合可以用一个五元组  来表示。式中各个参数含义如下：

1. N为HMM的隐状态个数。令符号 为t时刻的状态，则隐含状态序列可表示为  ，  。其中，  为隐状态集。
2. M为HMM的可观测状态个数。令符号 为t时刻的观测，则观测序列可表示为  ， 。其中，  为观测集。
3. A为N*N的状态转移概率分布矩阵。 ，表示状态由  转移到  的概率。
4. B为N*M的观测概率分布矩阵。 ，表示t时刻状态为  的条件下出现观察值  的概率。
5.  为初始状态分布。 ，表示序列初始状态为  的概率。

一个标准的HMM需要解决三个基本问题：

1. 评估问题，也被称为识别问题。给定模型  和观测序列O，计算每个HMM模型产生当前观测序列的输出概率  ，通过对比各个模型输出概率，概率最大的模型即为识别结果。
2. 解码问题。给定模型  和观测序列O，寻找最有可能产生观察序列O的隐含状态序列的过程，即寻找最佳状态序列 最好地解释观测序列O。
3. 学习问题。已知观测序列O，估计模型  参数，使得  最大。

二、问题描述：

三枚硬币，随机选择一枚，进行抛掷，记录抛掷结果。

该过程可以用隐马尔可夫模型来描述。假设得到一个HMM模型为，模型参数为

1. N=3，S={"coin1", "coin2", "coin3"}
2. M=2，V={H,T}。其中，H代表head，正面；T代表tail，反面
3. 
4. 
5. 

对应到HMM需要解决三个基本问题如下：

1.评估问题。求给定上述模型，观察到下列抛掷结果的概率是多少？

2.解码问题。求给定上述模型，若观察到上述抛掷结果，最可能的硬币选择序列（状态转换序列）是什么？

3.学习问题。若上述模型中的状态转移矩阵A，状态输出矩阵B和初始状态分布  均未知，如何根据问题1中的观察序列得到他们？

三、基于hmmlearn实现三个问题：

1.首先解决评估问题与解码问题，模型参数同问题描述。

建立模型

# -*- coding: utf-8 -*-
import numpy as np
from hmmlearn import hmm
import math

# 定义N
states = ["coin1", "coin2", "coin3"]
n_states = len(states)  # 获取状态个数

# 定义M
observations = ["Head", "Tail"]
n_observations = len(observations)  # 获取观测个数

# 定义pi
start_probability = np.array([0.3, 0.3, 0.4])

# 定义A
transition_probability = np.array([
  [0.5, 0.45, 0.05],
  [0.45, 0.1, 0.45],
  [0.45, 0.45, 0.1]
])

# 定义B
emission_probability = np.array([
  [0.5, 0.5],
  [0.75, 0.25],
  [0.25, 0.75]
])
# 建立HMM模型，共有n_states个隐状态
model = hmm.MultinomialHMM(n_components=n_states)
# HMM模型赋值
model.startprob_=start_probability
model.transmat_=transition_probability
model.emissionprob_=emission_probability

得到模型后，直接调用hmmlearn的内置函数即可求解。其中，输入的观测序列{0,0,1,0,1,1}对应  。

# 解码问题
seen = np.array([[0, 0, 1, 0, 1, 1]]).T
logprob, coin = model.decode(seen, algorithm="viterbi")
print coin

# 评估问题
probability = model.score(seen)
print math.e**probability

对于评估问题，模型返回的值是以自然对数为底的对数概率值，要得到真实概率值只需执行math.e**probability即可。

2.学习问题，即已知两个观测序列  ，N=3，M=2，求出A,B,  。

建立模型

# -*- coding: utf-8 -*-
import numpy as np
from hmmlearn import hmm
import math

states = ["coin1", "coin2", "coin3"]
n_states = len(states)

observations = ["Head", "Tail"]
n_observations = len(observations)
model = hmm.MultinomialHMM(n_components=n_states, n_iter=20, tol=0.01)
# n_inter: 要执行的最大迭代数
# tol: 收敛阈值。如果log-likelihood的增益低于这个值，EM就会停止。

输入观测序列，训练模型，并输出概率值

# 观测序列
X = np.array([[0, 0, 1, 0, 1, 1]])
# 模型训练
model.fit(X2)
print model.startprob_
print model.transmat_
print model.emissionprob_
print math.e**model.score(X)

多次训练模型，从中找到一组相对最优的参数值。