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次元目, ...) となっている. なんでこっちで１次元に落とさないといけないんだろう
. 内部的にやってくれればいいのにー.
まあこれで面倒な部分は大体書いた筈. 今回お世話になったページ.
