HMM求解三大问题的小例子

关于HMM求解三大问题的小例子，下面两个博客讲的不错:
https://www.jianshu.com/p/0f175b9781de
用hmmlearn学习隐马尔科夫模型HMM

笔记：
hmmlearn实现了三种HMM模型类，按照观测状态是连续状态还是离散状态，可以分为两类。GaussianHMM和GMMHMM是连续观测状态的HMM模型，而MultinomialHMM是离散观测状态的模型

对于连续观测状态的HMM模型，GaussianHMM类假设观测状态符合高斯分布，而GMMHMM类则假设观测状态符合混合高斯分布。一般情况下我们使用GaussianHMM即高斯分布的观测状态即可。

如果观测序列是一维的，则观测状态的概率密度函数是一维的普通高斯分布。如果观测序列是N维的，则隐藏状态对应的观测状态的概率密度函数是N维高斯分布。高斯分布的概率密度函数参数可以用μ表示高斯分布的期望向量，Σ表示高斯分布的协方差矩阵。在GaussianHMM类中，“means”用来表示各个隐藏状态对应的高斯分布期望向量μ形成的矩阵，而“covars”用来表示各个隐藏状态对应的高斯分布协方差矩阵Σ形成的三维张量

协方差矩阵的概念复习：
http://blog.sina.com.cn/s/blog_4aa4593d01012am3.html
https://blog.csdn.net/itplus/article/details/11452743

简单理解：一维单高斯 多维单高斯 混合多高斯GMM
https://blog.csdn.net/lzrtutu/article/details/77783081
详解EM算法与混合高斯模型(Gaussian mixture model, GMM)
https://blog.csdn.net/lin_limin/article/details/81048411

笔记：

对角矩阵(diagonal matrix)是一个主对角线之外的元素皆为0的矩阵

对称矩阵（Symmetric Matrices）是指元素以主对角线为对称轴对应相等的矩阵。在线性代数中，对称矩阵是一个方形矩阵，其转置矩阵和自身相等

服从二维高斯分布的数据主要集中在一个椭圆内部，服从三维的数据集中在一个椭球内部。

那么看过一维单高斯，多维单高斯，我们来看看经常被提到的GMM是个什么东东，GMM的中文意思其实就是高斯混合模型。单高斯意思是只有一个高斯，而GMM就是有多个高斯混合而成。字面意思可以理解有两个或者两个以上的高斯混合一起的可以被称为GMM。每个高斯实际就是一个component，而每个component可以是一个一维单高斯也可以是一个多维单高斯

所以，混合高斯模型并不是什么新奇的东西，它的本质就是融合几个单高斯模型，来使得模型更加复杂，从而产生更复杂的样本。理论上，如果某个混合高斯模型融合的高斯模型个数足够多，它们之间的权重设定得足够合理，这个混合模型可以拟合任意分布的样本。

下面主要整理一下代码：

1. 学习，也就是参数估计问题。已知观察序列，来对HMM模型的参数进行估计。

'''
下面抽取5次，得到已知的观察序列，
来对HMM的参数进行估计，即使用MultinomialHMM进行参数的训练

'''

import numpy as np
import hmmlearn.hmm as hmm

states = ['盒子1', '盒子2', '盒子3']
obs = ['白球', '黑球']
n_states = len(states)
m_obs = len(obs)

model2 = hmm.MultinomialHMM(n_components=n_states, n_iter=20, tol=0.001)
X2 = np.array([
    [0, 1, 0, 0, 1],
    [0, 0, 0, 1, 1],
    [1, 1, 0, 1, 0],
    [0, 1, 0, 1, 1],
    [0, 0, 0, 1, 0]
])
model2.fit(X2)
print("输出根据数据训练出来的π")
print(model2.startprob_)
print("输出根据数据训练出来的A")
print(model2.transmat_)
print("输出根据数据训练出来的B")
print(model2.emissionprob_)

输出结果如下：

输出根据数据训练出来的π
[  9.99999998e-01   1.81378454e-09   8.05262820e-30]
输出根据数据训练出来的A
[[ 0.30459785  0.18910843  0.50629372]
 [ 0.44920744  0.35079347  0.19999909]
 [ 0.48427612  0.45080508  0.0649188 ]]
输出根据数据训练出来的B
[[ 0.9600675   0.0399325 ]
 [ 0.48750091  0.51249909]
 [ 0.02464824  0.97535176]]

2. 解码问题。已知观察序列，求什么样的隐藏状态序列最可能生成一个给定的观察序列

'''
作者：白尒摩斯
链接：https://www.jianshu.com/p/0f175b9781de
来源：简书

'''

import numpy as np
import hmmlearn.hmm as hmm

# 首先定义变量
status = ['盒子1', '盒子2', '盒子3']   # 隐藏的状态集合
obs = ['白球', '黑球']    # 观察值集合
n_status = len(status)    # 隐藏状态的长度
m_obs = len(obs)          # 观察值的长度
start_probability = np.array([0.2, 0.5, 0.3])   
'''
初始概率分布： π 表示初次抽时，抽到1盒子的概率是0.2，抽到2盒子的概率是0.5，抽到3盒子的概率是0.3
'''
transition_probability = np.array([
    [0.5, 0.4, 0.1],
    [0.2, 0.2, 0.6],
    [0.2, 0.5, 0.3]
])
'''
状态转移概率矩阵 A[0][0]=0.5 表示当前我抽到1盒子，下次还抽到1盒子的概率是0.5
'''
emission_probalitity = np.array([
    [0.4, 0.6],
    [0.8, 0.2],
    [0.5, 0.5]
])
'''
观测概率矩阵 B：如最初的图，B[0][0]=0.4,表示第一个盒子抽到白球概率0.4，B[0][1]=0.6，表示第一个盒子抽到黑球概率0.6

以上各个参数解释请参考下面链接：
https://www.jianshu.com/p/da633461684f
https://www.jianshu.com/p/0f175b9781de
'''

# 下面开始定义模型
'''
hmmlearn中主要有两种模型，分布为：GaussianHMM和MultinomialHMM；
如果观测值是连续的，那么建议使用GaussianHMM，否则使用MultinomialHMM

参数：
初始的隐藏状态概率π参数为: startprob；
状态转移矩阵A参数为: transmat;
状态和观测值之间的转移矩阵B参数为: emissionprob_(MultinomialHMM模型中)或者在GaussianHMM模型中直接给定均值(means)和方差/协方差矩阵(covars)

'''

model = hmm.MultinomialHMM(n_components=n_status)  # 观测值，球的黑或白，是离散的，n_status隐藏状态的长度
model.startprob_ = start_probability
model.transmat_ = transition_probability
model.emissionprob_ = emission_probalitity


'''
下面运行viterbi预测问题。已知观察序列（白球 黑球 白球 白球 黑球），
求什么样的隐藏状态序列（盒子2 盒子3 盒子2 盒子2 盒子3）最可能生成一个给定的观察序列。

status = ['盒子1', '盒子2', '盒子3']
obs = ['白球', '黑球'] 

'''

se = np.array([[0, 1, 0, 0, 1]]).T   # （白球 黑球 白球 白球 黑球）
logprob, box_index = model.decode(se, algorithm='viterbi')
print("颜色:", end="")
print(" ".join(map(lambda t: obs[t], [0, 1, 0, 0, 1])))
print("盒子:", end="")
print(" ".join(map(lambda t: status[t], box_index)))
print("概率值:", end="")
print(np.exp(logprob)) # 这个是因为在hmmlearn底层将概率进行了对数化，防止出现乘积为0的情况

输出结果如下：

颜色:白球 黑球 白球 白球 黑球
盒子:盒子2 盒子3 盒子2 盒子2 盒子3
概率值:0.002304

3. 评估问题。这个暂时没有找到例子，但是看了上面两个例子，大概知道怎么做。待补充。