HMM隐马尔可夫模型来龙去脉（一）

目录

隐马尔可夫模型HMM学习导航

　　一、认识贝叶斯网络

　　　　1、概念原理介绍

　　　　2、举例解析

　　二、马尔可夫模型

　　　　1、概念原理介绍

　　　　2、举例解析

　　三、隐马尔可夫模型

　　　　1、概念原理介绍

　　　　2、举例解析

　　四、隐马尔可夫模型简单实现

　　五、完整代码

　　六、结语

隐马尔可夫模型HMM学习导航

NLP学习记录，这一章从概率图模型开始，学习常见的图模型具体的原理以及实现算法，包括了有向图模型：贝叶斯网络(BN)、（隐）马尔可夫模型(MM/HMM)，无向图模型：马尔可夫网络(MN)、条件随机场(CRF)。学习前提条件需要一定的概率论与数理统计知识，里面许多方法都是概率统计知识。

图模型框架：

一、认识贝叶斯网络

1、概念原理介绍

贝叶斯网络是一种基于概率推理的数学模型，属于有向无环图，理论基础就是我们熟悉的贝叶斯公式。

formula  

公式中，事件Bi的概率为P(Bi)，事件Bi已发生条件下事件A的概率为P(A│Bi)，事件A发生条件下事件Bi的概率为P(Bi│A)，公式分母就是全概率公式。

构造贝叶斯网络任务复杂，包括表示、推断、学习三个方面的问题。这里就简单介绍。

  1. 表示：某个随机变量集合\large x=\{X_1,... ,X_n\}上联合概率分布P，贝叶斯网络的表示就是在随机变量仅有两种取值的简单情况下，对 \large 2^n 种不同取值的联合概率分布进行说明，可以看出，这是计算量巨大的任务。
  2. 推断：贝叶斯网络是变量及其关系的完整模型，可以从观察到或已知的变量推断另一些变量的变化，这种根据已知的证据计算变量的后验分布称为概率推理。
  3. 学习：指的是BN中参数的学习，目的是决定变量之间的关联量化关系（依存强度估计）。参数学习方法包括：最大似然估计法、最大后验概率法、期望最大化法和贝叶斯估计法（使用最多）。

2、举例解析

从前面一片理论介绍和概念中，似乎很不直观，还没理解贝叶斯网络到底是什么，现在我们通过简单例子更加直观深刻地去认识。

假设一篇关于南海岛屿的新闻（News），新闻里可能包含关于介绍南海岛屿历史的内容（History），还有介绍旅游风光的内容（Sightseeing）。现在对此构造建立一个简单的贝叶斯网络。

现在用N、H、S表示这3个事件，每个事件存在两种可能，记为T（包含/是）、F（不包含/不是），通过假设概率建模：

New
T F
0.2 0.8
  Sightseeing
News T F
F 0.4 0.6
T 0.1 0.9
    History
Sightseeing News T F
F F 0.5 0.5
F T 0.3 0.7
T F 0.8 0.2
T T 0.4 0.6

一个简单的贝叶斯网络：

上面表格枚举了所有可能的情况概率分布，现在根据表格内容和建立的BN图关系，三个事件的联合概率：

即 \large P(H,S,N)=P(H|S,N)\times P(S|N)\times P(N) .

给出一个问题：如果一篇新闻中含有南海岛屿历史相关的内容，那么该新闻是关于南海新闻的可能性是多少呢？

显然，模型能够解决这个问题，根据题意我们得出，就是计算 H=T的情况下 N=T 的条件概率 \large P(N=T|H=T)。

进一步求解得到

\large P(N=T|H=T)=\frac{P(N=T,H=T)}{P(H=T)}= \frac{\sum_{S\in {T,F}}P(H=T,S,N=T)}{\sum_{N,S\in {T,F}}P(H=T,S,N)}

现在直接代入表格的数据，按H,S,N顺序，分子为TTT+TFT，分母TTT+TFT+TTF+TFF. 由此计算得到：

\large \frac{0.4\times 0.1\times 0.2+0.3\times 0.9\times 0.2}{0.4\times 0.1\times 0.2+0.3\times 0.9\times 0.2+0.8\times0.4\times0.8+0.5\times0.6\times0.8} =  \large 11.11\% .

二、马尔可夫模型

1、概念原理介绍

马尔可夫模型是一个可见的过程，隐马尔可夫模型则是将随机事件的状态转移过程屏蔽，所以学习马尔可夫模型作为前提。

马尔可夫模型（Markov Model）描述了一种随机函数，是随着时间不断变化的过程。换言之，我们观察到每个随机变量序列并不是相互独立的，而是依赖于其前面序列的状态。

2、举例解析

假设系统有N个状态S=\{s_1,s_2,...,s_N\}, 一个随机序列Q=\{q_1,q_2,...,q_N\}, t 时刻状态q_t=s_j，那它的概率如何计算呢？

计算t时刻状态的概率需要给出前面t-1个状态的关系，概率表示为：

P(q_t=s_j|q_{t-1}=s_i,q_{t-2}=s_k,......) ,相当于求解条件概率。

对于马尔可夫模型来说，有这样的特定条件，系统t时刻的状态只与其在时间t-1状态相关，所以表示为如下，称为一阶马尔可夫链。

P(q_t=s_j|q_{t-1}=s_i,q_{t-2}=s_k,......)=P(q_t=s_j|q_{t-1}=s_i).

进而由上面的是公式得到一个随机过程表示，构成状态转移矩阵。

P(q_t=s_j|q_{t-1}=s_i)=a_{ij}(1\leq i,j\leq N), 满足a_{ij}\geqslant 0,\sum_{j=1}^{N}a_{ij}=1 . 可见N个状态可以构成N*N个状态转移概率。

下面用更加直观例子解析马尔可夫模型。

假设一段文字中有三类词性：名词、动词、形容词。这三类词性的出现可以表示为马尔可夫模型的三个状态，描述为：

状态a:名词

状态b:动词

状态c:形容词

假设转移矩阵如下：

A=[a_{ij}]=\begin{bmatrix} 0.3 &0.5 &0.2 \\ 0.5 &0.3& 0.2 \\ 0.4 &0.2 &0.4 \end{bmatrix}

各状态关系，每个状态所有发射弧概率之和为1

根据此模型M，计算某个句子序列词性出现顺序 O = "名动形名"的概率：

P(O|M)=P(a,b,c,a|M)=P(a)\times P(b|a)\times P(c|b)\times P(a|c) = 1\times a_{12} \times a_{23} \times a_{31} =0.5 \times 0.2\times 0.4=0.04

三、隐马尔可夫模型

1、概念原理介绍

在前面马尔可夫模型，每个状态都是可见的。而隐马尔可夫模型，就是将这状态序列隐蔽，变为状态转移过程不可见。HMM是一个双重随机过程。它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。

2、举例解析

一个暗室里面有N个口袋，每个口袋有M个不同颜色小球。这个过程包括实验者和观察者。实验者每次根据某个概率分布随机选取一个口袋，再根据某个概率分布从这个口袋中取出一个小球，并向观察者展示该球的颜色。

我们理解一下这个过程，观察者只能知道每次取出球的颜色序列，而不知是从哪个口袋取出的，这个口袋序列的状态转移只有暗室里的实验者知道。所以，观察者的任务就是根据可观察序列（球颜色序列）来推断出隐藏的转移过程。

由此，总结出一个HMM有以下组成：

  1. 模型状态数目N（如例子中的口袋数量）
  2. 每个状态可能输出的不同符号数目M（一个口袋不同颜色球的数目）
  3. 状态转移矩阵A={a_{ij}} (a_{ij}是实验者从一口袋 i 转向口袋 j 的概率)。P(q_t=s_j|q_{t-1}=s_i)=a_{ij}(1\leq i,j\leq N)，a_{ij}\geqslant 0,\sum_{j=1}^{N}a_{ij}=1。
  4. 观察到符号的概率（符号发射概率）分布矩阵B={b_{j}(k)} (表示实验者从第j个袋取出第k种颜色的球)。
  5. 初始状态概率分布π。

四、隐马尔可夫模型简单实现

初始化状态数目n和样本数量，得到初始概率：

import numpy as np
from hmmlearn import hmm
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.metrics.pairwise import pairwise_distances_argmin
 
np.random.seed(28)
n = 5  # 隐状态数目
n_samples = 500  # 样本数量
pi = np.random.rand(n)
pi /= pi.sum()
print('初始概率：')
print(pi)

计算转移矩阵：

A = np.random.rand(n, n)
mask = np.zeros((n, n), dtype=np.bool)
mask[0][1] = mask[0][4] = True
mask[1][0] = mask[1][2] = True
mask[2][1] = mask[2][3] = True
mask[3][2] = mask[3][4] = True
mask[4][0] = mask[4][3] = True
A[mask] = 0
for i in range(n):
    A[i] /= A[i].sum()
print('转移概率:')
print(A)

给定均值和方差：

# 给定均值
means = np.array(((30, 30, 30), (0, 50, 20), (-25, 30, 10),
                  (-15, 0, 25), (15, 0, 40)), dtype=np.float)
for i in range(n):
    means[i, :] /= np.sqrt(np.sum(means ** 2, axis=1))[i]
covars = np.empty((n, 3, 3))
for i in range(n):
    covars[i] = np.diag(np.random.rand(3) * 0.02 + 0.001)

产生模拟数据、模型构建及估计参数：

# 产生对应的模拟数据
model = hmm.GaussianHMM(n_components=n, covariance_type='full')
model.startprob_ = pi
model.transmat_ = A
model.means_ = means
model.covars_ = covars
sample, labels = model.sample(n_samples=n_samples, random_state=0)

# 模型构建及估计参数
model = hmm.GaussianHMM(n_components=n, n_iter=10)
model.fit(sample)
y = model.predict(sample)
np.set_printoptions(suppress=True)
# print('##估计初始概率：')
# print(model.startprob_)
# print('##估计转移概率：')
# print(model.transmat_)
# print('##估计均值：\n')
# print(model.means_)
# print('##估计方差：\n')
# print(model.covars_)

根据类别信息更改顺序：

# 根据类别信息更改顺序
order = pairwise_distances_argmin(means, model.means_, metric='euclidean')
# print(order)
pi_hat = model.startprob_[order]
A_hat = model.transmat_[order]
A_hat = A_hat[:, order]
means_hat = model.means_[order]
covars_hat = model.covars_[order]
change = np.empty((n, n_samples), dtype=np.bool)
for i in range(n):
    change[i] = y == order[i]
for i in range(n):
    y[change[i]] = i
acc = np.mean(labels == y) * 100
print('准确率：%.2f%%' % acc)

3D画图分类结果：

mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False
fig = plt.figure(figsize=(8, 8), facecolor='w')
ax = fig.add_subplot(111, projection='3d')  #
colors = plt.cm.Spectral(np.linspace(0, 1, n))
ax.scatter(sample[:, 0], sample[:, 1], sample[:, 2], s=50, c=labels,
           cmap=plt.cm.Spectral, marker='o', label=u'观测值', depthshade=True)
plt.plot(sample[:, 0], sample[:, 1], sample[:, 2], lw=0.1, color='#A07070')
colors = plt.cm.Spectral(np.linspace(0, 1, n))
ax.scatter(means[:, 0], means[:, 1], means[:, 2], s=300, c=colors,
           edgecolor='r', linewidths=1, marker='*', label=u'中心')
 
x_min, y_min, z_min = sample.min(axis=0)
x_max, y_max, z_max = sample.max(axis=0)
x_min, x_max = expand(x_min, x_max)
y_min, y_max = expand(y_min, y_max)
z_min, z_max = expand(z_min, z_max)
ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))
ax.set_zlim((z_min, z_max))
plt.legend(loc='upper left')
plt.grid(True)
plt.title(u'GMHMM参数估计和类别判定', fontsize=12)
plt.show()
View Code

五、完整代码

import numpy as np
from hmmlearn import hmm
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from sklearn.metrics.pairwise import pairwise_distances_argmin
 
 
def expand(a, b):
    return 1.05 * a - 0.05 * b, 1.05 * b - 0.05 * a
 
 
np.random.seed(28)
n = 5  # 隐状态数目
n_samples = 500  # 样本数量
pi = np.random.rand(n)
pi /= pi.sum()
# print('初始概率：')
# print(pi)
 
A = np.random.rand(n, n)
mask = np.zeros((n, n), dtype=np.bool)
mask[0][1] = mask[0][4] = True
mask[1][0] = mask[1][2] = True
mask[2][1] = mask[2][3] = True
mask[3][2] = mask[3][4] = True
mask[4][0] = mask[4][3] = True
A[mask] = 0
for i in range(n):
    A[i] /= A[i].sum()
# print('转移概率:')
# print(A)
 
# 给定均值
means = np.array(((30, 30, 30), (0, 50, 20), (-25, 30, 10),
                  (-15, 0, 25), (15, 0, 40)), dtype=np.float)
for i in range(n):
    means[i, :] /= np.sqrt(np.sum(means ** 2, axis=1))[i]
# print(means)
 
# 给定方差
covars = np.empty((n, 3, 3))
for i in range(n):
    covars[i] = np.diag(np.random.rand(3) * 0.02 + 0.001)  # np.random.rand ∈[0,1)
# print(covars)
 
# 产生对应的模拟数据
model = hmm.GaussianHMM(n_components=n, covariance_type='full')
model.startprob_ = pi
model.transmat_ = A
model.means_ = means
model.covars_ = covars
sample, labels = model.sample(n_samples=n_samples, random_state=0)
 
# 模型构建及估计参数
model = hmm.GaussianHMM(n_components=n, n_iter=10)
model.fit(sample)
y = model.predict(sample)
np.set_printoptions(suppress=True)
 
# 根据类别信息更改顺序
order = pairwise_distances_argmin(means, model.means_, metric='euclidean')
# print(order)
pi_hat = model.startprob_[order]
A_hat = model.transmat_[order]
A_hat = A_hat[:, order]
means_hat = model.means_[order]
covars_hat = model.covars_[order]
change = np.empty((n, n_samples), dtype=np.bool)
for i in range(n):
    change[i] = y == order[i]
for i in range(n):
    y[change[i]] = i
acc = np.mean(labels == y) * 100
print('准确率：%.2f%%' % acc)
 
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False
fig = plt.figure(figsize=(8, 8), facecolor='w')
ax = fig.add_subplot(111, projection='3d')  #
colors = plt.cm.Spectral(np.linspace(0, 1, n))
ax.scatter(sample[:, 0], sample[:, 1], sample[:, 2], s=50, c=labels,
           cmap=plt.cm.Spectral, marker='o', label=u'观测值', depthshade=True)
plt.plot(sample[:, 0], sample[:, 1], sample[:, 2], lw=0.1, color='#A07070')
colors = plt.cm.Spectral(np.linspace(0, 1, n))
ax.scatter(means[:, 0], means[:, 1], means[:, 2], s=300, c=colors,
           edgecolor='r', linewidths=1, marker='*', label=u'中心')
 
x_min, y_min, z_min = sample.min(axis=0)
x_max, y_max, z_max = sample.max(axis=0)
x_min, x_max = expand(x_min, x_max)
y_min, y_max = expand(y_min, y_max)
z_min, z_max = expand(z_min, z_max)
ax.set_xlim((x_min, x_max))
ax.set_ylim((y_min, y_max))
ax.set_zlim((z_min, z_max))
plt.legend(loc='upper left')
plt.grid(True)
plt.title(u'GMHMM参数估计和类别判定', fontsize=12)
plt.show()
View Code

六、结语

通过这篇学习记录，我们初步认识了HMM隐马尔可夫模型的具体内部逻辑，另外举例解析和最后python简单实现一个例子，更直观的去理解什么是HMM，它是如何工作的。到这里可能觉得代码实现还有部分陌生的方法，因为HMM是一个比较复杂的任务，除了本文简单入门，HMM还需要解决三个基本问题，认真阅读的朋友会看到在代码中有所体现：估计问题、序列问题和训练问题（参数估计）。

这些问题等待接下来学习之后再来具体介绍，期待下一篇更核心的内容解析。

参考资料：《统计自然语言处理》、https//www.jianshu.com/p/083c8dfb9f0a

我的博客园：https://www.cnblogs.com/chenzhenhong/p/13537680.html

我的CSDN: https://blog.csdn.net/Charzous/article/details/108111860

 

 

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Charzous/article/details/108111860

 

posted @ 2020-08-20 21:44  Charzueus  阅读(1273)  评论(0)  编辑  收藏  举报