PythonからのHMMライブラリの利用

PythonからHMM(隠れマルコフモデル)ライブラリであるghmmを利用する方法の紹介です.
インストール方法は公式ページhttp://ghmm.org/を参照.

次のサンプルコードは, いかさま師をモチーフにしたHMMを生成しています. HMMの出力はいかさま師のサイコロの目に相当します. イカサマ師は2種類のサイコロ(普通の奴と6しか出ない奴)を適当に切り替えて使っています. その出力列だけを見て, 状態遷移確率と出力分布の確率を学習しているのが以下のコードです.

import ghmm
import random

def random_matrix(n,m):
	"n,mのrandom配列を返す"
	result = [[0.0] * m for i in range(n)]

	for i in range(n):
		for j in range(m):
			result[i][j] = random.random()

	return result

#出力確率変数の値域
values = ghmm.IntegerRange(1,7)

#遷移行列
A = [[0.9,0.1],[0.3,0.7]]
efair = [1.0 / 6] * 6
eloaded = [0,0,0,0,0,1]
#状態ごとの出力分布
B = [efair,eloaded]
pi = [1.0,0.0]

m = ghmm.HMMFromMatrices(values,ghmm.DiscreteDistribution(values),A,B,pi)

#HMMからサンプリング(長さ1000のものを100個)
train_seqs = m.sample(100,1000)

print "correct:"
print m
print sum(m.viterbi(train_seqs)[1])

best_score = -10000000
best_HMM = None

for loop in range(10):
	initA = random_matrix(2,2)
	initB = random_matrix(2,6)

	detecter = ghmm.HMMFromMatrices(
		values,ghmm.DiscreteDistribution(values),
		initA,
		initB,
		random_matrix(1,2)[0])

	#遷移確率や出力分布が合計1になるように正規化
	detecter.normalize()
	#学習
	detecter.baumWelch(train_seqs)

	if best_score < sum(detecter.viterbi(train_seqs)[1]):
		best_score = sum(detecter.viterbi(train_seqs)[1])
		best_HMM = detecter

print best_HMM
print best_score

出力は次のようになり, だいたい正しくそれぞれの確率を推定できていることがわかります.

correct:
DiscreteEmissionHMM(N=2, M=6)
  state 0 (initial=1.00)
    Emissions: 0.17, 0.17, 0.17, 0.17, 0.17, 0.17
    Transitions: ->0 (0.90), ->1 (0.10)
  state 1 (initial=0.00)
    Emissions: 0.00, 0.00, 0.00, 0.00, 0.00, 1.00
    Transitions: ->0 (0.30), ->1 (0.70)

-166597.17182
DiscreteEmissionHMM(N=2, M=6)
  state 0 (initial=0.00)
    Emissions: 0.01, 0.01, 0.01, 0.01, 0.01, 0.96
    Transitions: ->0 (0.72), ->1 (0.28)
  state 1 (initial=1.00)
    Emissions: 0.17, 0.17, 0.17, 0.17, 0.17, 0.15
    Transitions: ->0 (0.11), ->1 (0.89)

-167385.703669

各関数の簡単な説明と, 注意すべき点を適当に. HMMFromMatricesは確率変数のドメインと, 状態遷移確率, 出力分布, 初期分布からHMMを生成します. この例だと出力分布に離散確率分布を扱っていますが,

#状態遷移確率
A = [[0.9,0.1],[0.3,0.7]]
#出力分布 (状態1の平均ベクトル, 状態1の共分散行列を1次元にしたもの), (状態2) ...
B = [[ [8.0,2.0] ,[1.0,0.0, 0.0,1.0] ],
	 [ [-3.0,-6.0] ,[1.0,0.0, 0.0,1.0] ] ]
pi = [1.0,0.0]
m = ghmm.HMMFromMatrices(values,ghmm.MultivariateGaussianDistribution(values),A,B,pi)

のようにしてやれば, 出力分布にガウス分布を用いることも出来ます. 混合ガウス分布も可能です.
ここでは共分散行列を1次元にしているのがポイントで, このライブラリはやたら2次元のものを1次元にさせたがります.

sample関数はHMMからのサンプリングです. viterbiは与えられた観測列に対する(確率が最も高い状態列, 尤度)を返します.
normalizeはHMMFromMatricesに与えた状態遷移確率などが確率分布として不正な場合(合計が1でない場合), 合計が1になるように正規化してくれる関数です.

baumWelchはBaumWelch法でHMMのパラメーターを調整してくれます. 注意すべき点は, 学習前のHMMの出力分布をそれぞれの状態で等しくするとゴミのような学習結果が出ます. 具体的には, 学習後もHMMの出力分布は各状態で等しくなります. このようにHMMは初期状態に悲しいぐらい依存します. 悲しい.

まあこんな感じで便利なライブラリです. 自分でPythonで実装しましたが, これに比べると悲しいぐらい遅いです.