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で実装しましたが, これに比べると悲しいぐらい遅いです.

ghmmにおける多次元連続確率分布の取り扱い

ghmmでろくに説明のない連続確率分布の取り扱いに関する解説. 大まかな説明は前回のページPythonからのHMMライブラリの利用 - Risky Duneを見ていただきたい.

以下のサンプルでは, 状態数2でそれぞれの状態において2次元の正規分布から値を出力するHMMを生成し, その観測列だけを用いてHMMのパラメーターを学習できるか実験している.

import ghmm
import random
import pylab
from scipy import *
import scipy.sparse as sparse
import scipy.linalg as linalg

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

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

	if n == 1:
		return result[0]
	return result

values = ghmm.Float()

A = [[0.9,0.1],[0.3,0.7]]
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)

train_seqs = m.sample(10,100)

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

best_score = -10000000
best_HMM = None

for loop in range(10):
	mu1 = random_matrix(2)
	mu2 = random_matrix(2)

	#分散共分散行列は半正定値
	sigma1 = random_matrix(2,2)
	L = sparse.tril(sigma1)
	#固有値が小さすぎるとバグるので
	sigma1 = dot(L,L.T).todense() * 100

	sigma2 = random_matrix(2,2)
	L = sparse.tril(sigma2)
	sigma2 = dot(L,L.T).todense() * 100

	sigma1 = sigma1[0].tolist()[0] + sigma1[1].tolist()[0]
	sigma2 = sigma2[0].tolist()[0] + sigma2[1].tolist()[0]
	
	initA = random_matrix(2,2)
	initB = [ [mu1,sigma1],[mu2,sigma2] ]

	detecter = ghmm.HMMFromMatrices(
		values,ghmm.MultivariateGaussianDistribution(values),
		initA,
		initB,
		random_matrix(2))

	detecter.normalize()
	detecter.baumWelch(train_seqs)

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

print best_HMM.getEmission(0,0)
print best_HMM.getEmission(1,0)
print best_score

このプログラムの出力は

correct:
([8.0, 2.0], [1.0, 0.0, 0.0, 1.0])
([-3.0, -6.0], [1.0, 0.0, 0.0, 1.0])
-3232.65593556
training:
([-2.9892426181496994, -5.9983553568542893], [0.97064467648541697, 0.022944396140138268, 0.022944396140138271, 1.0151665909275061])
([8.0205274775166817, 2.0527254408839655], [0.99776084732721204, -0.020820373620631539, -0.020820373620631539, 1.0776094241497327])
-3229.50292233

となり, だいたい平均と共分散行列は正しく推定できていることがわかる. では, 勘所を説明していこう.

以下は, HMMの初期共分散行列を生成している.

#分散共分散行列は半正定値
sigma1 = random_matrix(2,2)
L = sparse.tril(sigma1)
#固有値が小さすぎるとバグるので
sigma1 = dot(L,L.T).todense() * 100
#１次元に直す
sigma1 = sigma1[0].tolist()[0] + sigma1[1].tolist()[0]

うっかり初期値に正定値でない共分散行列を渡すと, 色んなパラメーターがnanになり動かない. そこで, 下三角行列を生成し, それを自身の転置と掛けることで半正定値行列を生成している.
さらに, 固有値が小さい行列を渡すと

GHMM ghmm.py:148 - sviterbi.c:ghmm_cmodel_viterbi(215): sequence can't be build from model

と連呼されるため, 100倍している(maxとか使った方がよかった).
そして, 共分散行列を1次元の形でうけとるため, 1次元に直している.

追記(10/04):
これは正確には固有値の問題では無い. 学習の過程で, 渡された観測列が現在のパラメーターでは, 生成できない際に

GHMM ghmm.py:148 - sviterbi.c:ghmm_cmodel_viterbi(215): sequence can't be build from model

というエラーが生じる. 例えば分散が0の共分散行列を持ってしまうとこれが生じるのはわかるだろう. よって学習のスタート時からある程度大きな分散を持つ必要がある. (k-meansアルゴリズムを用いて初期化するのが望ましい)

他は前回のページで述べた通りである.

最後にHMMを学習させる際に, HMMのサンプルからでは無く, 自分でリストから渡したい場合について説明しよう. HMMに観測列を渡す場合は, 1回の観測なら(たぶん)EmissionSequence, 複数の観測ならSequenceSetを用いる. 具体的にはこんな感じ

observed = [[0.0,4.0,0.1,4.1,-0.2,4.0,1.0,0.0,0.8,0.2],[0.0,4.0,0.1,4.1,-0.2,4.0,0.9,0.1]]
train_seqs = ghmm.SequenceSet(values,observed)
detecter.baumWelch(train_seqs)

リストの中身はobserved = (1番目の観測列), (2番目の観測列)であり, 更に(1番目の観測列) = (1番目の1回目に観測した値の1次元目, 1番目の1回目に観測した値の2次元目, 2回目に観測した値の1次元目, 2回目に観測した値の2次元目, ...) となっている. なんでこっちで１次元に落とさないといけないんだろう
. 内部的にやってくれればいいのにー.

まあこれで面倒な部分は大体書いた筈. 今回お世話になったページ.

• サンプルコードhttp://nullege.com/codes/search/ghmm.HMMFromMatrices.baumWelch
• リファレンスhttp://ghmm.sourceforge.net/Pydoc/ghmm.html

初めてのGreaseMonkey

Twitterのホームに, The Interviewsへのリンクを追加するGreasemonkeyを作成(http://userscripts.org/scripts/show/114020). Twitterのscreen_nameとThe InterviewsのIDが同じ場合, 名前の下にリンクが出現する.

初めてGreasemonkeyを書いたけど, 結構楽しい.

ちょっと混乱したポイントとしては, Twitterをブラウザで開いた状態で, URLバーにtwitter.com/user_nameみたいに入力して移動した場合, 新しくページが開かれるのではなく, ページの差分だけ読み込むため, Greasemonkeyが新規に読み込まれなかった所.
これは, もう3秒ごとにスクリプトを実行することで無理やり解決. まあそんなコストでも無いだろう.

もう一つ重要なのはXMLHTTPRequestはクロスドメインで呼び出せない点.
The Interviewsにアカウントが存在するか確認するためにはクロスドメインアクセスが必要なため, GM_xmlhttpRequestを用いることで解決した.
この関数はGoogle Chromeに存在しないので, 今はfirefoxでしか動作しない. 一応, The Interviewsへのアクセスをプロキシを介して行うことで, chromeでもクロスドメインでXMLHTTPRequestを動かせるらしいんだけど, 今回は見送った.

最後に, スクリプトを公開したんだけど, これが凄く楽だった. userscripts.orgさん, ぱねぇ.

参考にしたサイト

• greasemonkey入門メモ | tokenの日記 | スラド
• Firebugで作るGreasemonkeyスクリプト〜入門と実践(From Kanasan.JS) | Blog.37to.net
• http://d.hatena.ne.jp/alexam/20110810/1312975246