HMM
HMMの概要
例えば、K人の話者がいて、収録した音声データから誰が話しているのかを推定する問題を考えましょう。簡単のため、K人の話者は同時には発話しないものとし、収録した音声から音の到来方向に関係する情報が観測できるとします。図１は、K=３の場合の観測値の例で、３人の話者に対応する３つの分布の山が見られます。
時刻1からNまでの音の到来方向に関係する情報（特徴量）を観測値X={x1,...,xN}とし、これに対応した話者の状態（誰が話しているのか）を Z={z1,...,zN} で表すとします。時刻 n の状態ベクトルzn=[zn,1,...,,zn,K]の要素zn,k は、k番目の話者が発話している時に1、発話していない時に0の値をとります。zn は、隠れ状態(hidden state)あるいは、潜在変数(latent variable)と呼ばれます。
モデルの学習
XとZの同時確率密度は式(1)のようになります。ここで、p(z0|π) は、時刻0での初期状態を表す確率密度です。π は、データが観測される前に、どの状態にあるかを示す事前確率を表します。p(xn|zn,φ) は、尤度と呼ばれます。そのパラメータφは、各状態における確率分布のパラメータを表し、ガウス分布の場合は、式(2)のような平均値ベクトル μkと共分散行列 Σk で表されます。図１の例では、分布の３つの山を、それぞれガウス分布で近似し、その平均値と共分散が、(μk, Σk )であるということになります。
p(zn|zn-1,A)は、状態遷移確率と呼ばれます。このパラメータ A は、状態遷移行列と呼ばれ、その各要素は、図１に示すように、各状態間の遷移確率を表します。式(1)に登場した全てのパラメタをまとめて、モデルのパラメータと呼び、式(3)のように表します。
最尤法では、式(1)において、Zについて周辺分布をとった次式の観測値全体に対する尤度を最大化するよう、パラメータ θ を決定します。しかし、式(4)の右辺を陽に計算することは一般に困難なため、EMアルゴリズムを用いて、反復的に最尤法の最適解を求めます。
隠れ状態の推定
モデルのパラメータ θ が推定できたら、ビタビ・アルゴリズム(Viterbi algorithm)などを用いて、隠れ状態 Z={z1,...,zN} を推定します。
参考文献
C. M. Bishop, "Pattern recognition and machine learning," Springer, 2006
インストール
ドキュメント: https://hmmlearn.readthedocs.io/en/latest/
インストール：https://github.com/hmmlearn/hmmlearn
図１：観測値のヒストグラムの例
図２：HMMにおける状態遷移
図３：生成された状態遷移
図４：生成された分布
図５：推定された状態遷移
図６：推定された分布
サンプルプログラム
hmmlearnには、
(a)モデルパラメータを与えて観測値を生成する、
(b)観測値からモデルパラメータを学習する、
(c)学習したモデルパラメータと観測値から隠れ状態を推定する、
の３つの機能があります。ここでは、この３つを使ってみましょう。まず、初期確率π、ガウス分布の平均値ベクトルと共分散行列{μ,Σ}、遷移行列A を与えて、観測値 X と 隠れ状態 Z を生成します。図３は生成された状態遷移、図４は、モデルパラメータに基づいた２次元ガウス分布とサンプルの散布図を表します。続いてメソッド fit() を用いて、パラメタの推定を行います。図６は推定されたガウス分布のモデルパラメタを用いて得られたガウス分布の形状を表しています。最後にメソッド predict() で隠れ状態の推定を行います。図５は、推定された隠れ状態を表しています。なお、見やすさのため、推定されたモデルパラメータと生成の時に用いたモデルパラメタを比較し、推定された隠れ状態の交換を行なっています（実際の応用で、この操作を行うには、先見的な知識が必要であり、HMMによる学習だけで交換の任意性を解決することはできません）。
import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt
from hmmlearn import hmm
# 観測値の生成
K = 3;
Nsample = 100
model = hmm.GaussianHMM(n_components=K, covariance_type="full")
model.startprob_ = np.array([0.3,0.5,0.2])
model.means_ = np.array([[-0.5,1.5], [1.5,0.3], [-0.8,-0.8]])
model.covars_ = np.array([[[0.3,-0.2],[-0.2,0.3]],
[[0.3, 0.0],[ 0.0,0.3]],
[[0.4, 0.2],[ 0.2,0.4]]])
model.transmat_ = np.array([[0.8,0.1,0.1],
[0.1,0.8,0.1],
[0.1,0.1,0.8]])
X,Z = model.sample(Nsample)
# 観測値のグラフ
Colors = ('b','g','r')
gxmin,gxmax = (-3,3)
gymin,gymax = (-3,3)
fig = plt.figure(figsize=(6,6))
for k,(mean,cov,color) in enumerate(zip(
model.means_,model.covars_,Colors)):
# 散布図
plt.scatter(X[Z == k,0], X[Z == k,1], color=color)
# 等高線
x, y = np.mgrid[gxmin:gxmax:0.01, gymin:gymax:0.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x;
pos[:, :, 1] = y
rv = st.multivariate_normal(mean,cov)
P = rv.pdf(pos)
plt.contour(x, y, P, 5)
plt.xlim(gxmin,gxmax)
plt.ylim(gymin,gymax)
plt.xlabel('feature 1')
plt.ylabel('feature 2')
plt.title('Generated')
# 状態遷移のグラフ
z = Z+1 # 生成された状態は{0,1,2}なので１を加える
t = np.arange(Nsample)
plt.figure(figsize=(6,4))
plt.plot(t,z)
for k in range(K):
plt.scatter(t[Z == k], z[Z == k], color=Colors[k])
plt.yticks([1,2,3])
plt.xlabel('n')
plt.ylabel('state')
plt.title('Generated')
# HMMパラメタの学習
remodel = hmm.GaussianHMM(n_components=K, covariance_type="full", n_iter=100)
remodel.fit(X)
print('[HMM] Mean:' )
print(remodel.means_)
print('[HMM] Covariance:')
print(remodel.covars_)
print('[HMM] Transition:')
print(remodel.transmat_)
# 隠れ状態の推定
Zestimate = remodel.predict(X)
# 状態の交換
permutation = np.zeros(K,dtype=int)
for k in range(K):
permutation[k] = int(np.argmin(np.linalg.norm(np.abs(model.means_-remodel.means_[k]),axis=1)))
# 推定結果のグラフ
Colors = ('b','g','r')
gxmin,gxmax = (-3,3)
gymin,gymax = (-3,3)
fig = plt.figure(figsize=(6,6))
for k,(mean,cov,color) in enumerate(zip(
remodel.means_,remodel.covars_,Colors)):
kk = permutation[k]
# 散布図
plt.scatter(X[Zestimate == k,0], X[Zestimate == k,1], color=Colors[kk])
# 等高線
x, y = np.mgrid[gxmin:gxmax:0.01, gymin:gymax:0.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x;
pos[:, :, 1] = y
rv = st.multivariate_normal(mean,cov)
P = rv.pdf(pos)
plt.contour(x, y, P, 5)
plt.xlim(gxmin,gxmax)
plt.ylim(gymin,gymax)
plt.xlabel('feature 1')
plt.ylabel('feature 2')
plt.title('Estimated')
# 推定された状態遷移のグラフ
z = permutation[Zestimate]+1
t = np.arange(Nsample)
plt.figure(figsize=(6,4))
plt.plot(t,z)
for k in range(K):
kk = permutation[k]
plt.scatter(t[Zestimate == k], z[Zestimate == k], color=Colors[kk])
plt.yticks([1,2,3])
plt.xlabel('n')
plt.ylabel('state')
plt.title('Estimated')
plt.show()