HMM隱馬爾可夫模型前因後果（一）

表示：某個隨機變量集合 $\large x=\{X_1,... ,X_n\}$ 上聯合機率分佈P，貝葉斯網絡的表示就是在隨機變量僅有兩種取值的簡單狀況下，對 $\large 2^n$ 種不一樣取值的聯合機率分佈進行說明，能夠看出，這是計算量巨大的任務。
推斷：貝葉斯網絡是變量及其關係的完整模型，能夠從觀察到或已知的變量推斷另外一些變量的變化，這種根據已知的證據計算變量的後驗分佈稱爲機率推理。
學習：指的是BN中參數的學習，目的是決定變量之間的關聯量化關係（依存強度估計）。參數學習方法包括：最大似然估計法、最大後驗機率法、指望最大化法和貝葉斯估計法（使用最多）。

二、舉例解析

從前面一片理論介紹和概念中，彷佛很不直觀，還沒理解貝葉斯網絡究竟是什麼，如今咱們經過簡單例子更加直觀深入地去認識。

假設一篇關於南海島嶼的新聞（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\%$ .

2、馬爾可夫模型

一、概念原理介紹

馬爾可夫模型是一個可見的過程，隱馬爾可夫模型則是將隨機事件的狀態轉移過程屏蔽，因此學習馬爾可夫模型做爲前提。

馬爾可夫模型（Markov Model）描述了一種隨機函數，是隨着時間不斷變化的過程。換言之，咱們觀察到每一個隨機變量序列並非相互獨立的，而是依賴於其前面序列的狀態。

二、舉例解析

假設系統有N個狀態 $S=\{s_1,s_2,...,s_N\}$ , 一個隨機序列 $Q=\{q_1,q_2,...,q_N\}$ , t 時刻狀態，那它的機率如何計算呢？

計算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$

3、隱馬爾可夫模型

一、概念原理介紹

在前面馬爾可夫模型，每一個狀態都是可見的。而隱馬爾可夫模型，就是將這狀態序列隱蔽，變爲狀態轉移過程不可見。HMM是一個雙重隨機過程。它用來描述一個含有隱含未知參數的馬爾可夫過程。其難點是從可觀察的參數中肯定該過程的隱含參數。

二、舉例解析

一個暗室裏面有N個口袋，每一個口袋有M個不一樣顏色小球。這個過程包括實驗者和觀察者。實驗者每次根據某個機率分佈隨機選取一個口袋，再根據某個機率分佈從這個口袋中取出一個小球，並向觀察者展現該球的顏色。

咱們理解一下這個過程，觀察者只能知道每次取出球的顏色序列，而不知是從哪一個口袋取出的，這個口袋序列的狀態轉移只有暗室裏的實驗者知道。因此，觀察者的任務就是根據可觀察序列（球顏色序列）來推斷出隱藏的轉移過程。

由此，總結出一個HMM有如下組成：

模型狀態數目N（如例子中的口袋數量）
每一個狀態可能輸出的不一樣符號數目M（一個口袋不一樣顏色球的數目）
狀態轉移矩陣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$ 。
觀察到符號的機率（符號發射機率）分佈矩陣B={ $b_{j}(k)$ } (表示實驗者從第j個袋取出第k種顏色的球)。
初始狀態機率分佈π。

4、隱馬爾可夫模型簡單實現

初始化狀態數目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

5、完整代碼

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

6、結語

經過這篇學習記錄，咱們初步認識了HMM隱馬爾可夫模型的具體內部邏輯，另外舉例解析和最後python簡單實現一個例子，更直觀的去理解什麼是HMM，它是如何工做的。到這裏可能以爲代碼實現還有部分陌生的方法，由於HMM是一個比較複雜的任務，除了本文簡單入門，HMM還須要解決三個基本問題，認真閱讀的朋友會看到在代碼中有所體現：估計問題、序列問題和訓練問題（參數估計）。

這些問題等待接下來學習以後再來具體介紹，期待下一篇更核心的內容解析。

參考資料：《統計天然語言處理》、https//www.jianshu.com/p/083c8dfb9f0a

個人博客園：http://www.javashuo.com/article/p-wentklzh-mo.html

個人CSDN: http://www.javashuo.com/article/p-egldhizj-nt.html

本文連接：http://www.javashuo.com/article/p-egldhizj-nt.html