Python 隠れマルコフモデル用ライブラリ hmmlearn の使い方メモ

隠れマルコフモデル (HMM; Hidden Markov Model) を実装した Python のライブラリ hmmlearn の使い方を理解したのでメモしておく。

HMM で扱う問題は3種類あって、それを理解していないと「使ってみたけどよくわからない」状態になりかねないので、まずはそれらをおさらいして、その後にそれぞれの問題に対応させながら API の使い方を説明していく。

HMM が扱う問題の種類

HMM では次の3つの種類の問題を扱う。

  1. 評価: 既知のパラメータで構成された HMM の出力として観測系列 \boldsymbol{x} が得られる確率 (尤度) を評価する
  2. 復号: 既知のパラメータで構成された HMM の出力として観測系列 \boldsymbol{x} が得られる確率が最も高くなる状態系列 \boldsymbol{s} を復号する
  3. 推定: 未知のパラメータの HMM からの出力として観測系列 \boldsymbol{x} が得られたとき、未知のパラメータを推定する

インストール

$ pip install hmmlearn
$ python 
>>> from hmmlearn import hmm
>>> hmm.GaussianHMM
<class 'hmmlearn.hmm.GaussianHMM'>

モデルを構築する

HMM には次の3つのパラメータがある。

  • 初期状態確率 \boldsymbol{\pi}
  • 遷移確率 \boldsymbol{A}
  • 出力確率 \boldsymbol{B}

これらのパラメータは学習で推定して指定することもできるが、まずは明示的に指定してモデルを構築する。

# ライブラリのインポート
import numpy as np
from hmmlearn import hmm

# GaussianHMM クラスから出力確率が正規分布に従う隠れマルコフモデルを作る。
# n_components パラメータで隠れ状態が3つあることを指定している。
model = hmm.GaussianHMM(n_components=3, covariance_type="full")

# 初期状態確率 π を指定する。
model.startprob_ = np.array([0.6, 0.3, 0.1])

# 遷移確率 A を指定する。
model.transmat_ = np.array([
        [0.7, 0.2, 0.1],
        [0.3, 0.5, 0.2],
        [0.3, 0.3, 0.4]])

# 出力確率 B を指定する。
# ただし出力は正規分布に従うと仮定しているため、正規分布のパラメータの
# 平均値 μ (means_) と、共分散 σ^2 (covars_) を指定する。
model.means_ = np.array([
        [0.0, 0.0],
        [10.0, 10.0],
        [100.0, 100.0]])
model.covars_ = 2 * np.tile(np.identity(2), (3, 1, 1))

ここでちょっとわかりにくいのは出力関数の正規分布のパラメータとして指定する model.means_model.covars_ だと思う。上の例においては、model.means_ によって

  • 状態 s_0 からは [0.0, 0.0]
  • 状態 s_1 からは [10.0, 10.0]
  • 状態 s_2 からは [100.0, 100.0]

を平均値とし、さらに model.covars_ に指定した共分散行列に従う正規乱数を出力することを定義している。

出力信号に2次元ベクトル (配列) を指定しているので、共分散行列 model.covars_ に指定する値には np.identity(2)ℝ^{2 \times 2} の行列を作っている。ここのサイズが一致しないとエラーになる。また1次元の場合であっても [[0], [10], [100]] のように配列で指定しないとエラーになる。

観測系列のサンプルを出力する

上で構築したモデルから sample メソッドを用いて観測系列サンプルを出力する。

# サンプル信号を 10 個出力する。
X, Z = model.sample(10)

戻り値は (観測系列, 状態系列) のタプルで得られる。その中身はこんな感じ。

>>> X
array([[  0.75706838,  -0.1280334 ],
       [ 10.3137587 ,  10.59635189],
       [ 99.27435882,  99.63294557],
       [ 98.74324108,  99.25505532],
       [ -0.5194227 ,   2.59958429],
       [  2.66607548,  -0.41898544],
       [  1.49013197,  -1.51010137],
       [  9.394274  ,  11.27156992],
       [ 10.1477135 ,   9.76125003],
       [  9.88244708,   9.33941603]])
>>> Z
array([0, 1, 2, 2, 0, 0, 0, 1, 1, 1])

観測値と状態の組み合わせを確認して、それらしい値が出力されていることが分かる。

[評価] 既知のモデルから観測系列 \boldsymbol{x} が得られる尤度を評価する

構築したモデルの score メソッドに観測系列 \boldsymbol{x} を渡すことで、観測系列の対数尤度が得られる。

上で出力した観測系列サンプルの対数尤度を評価すると

>>> model.score(X)
-41.160325787201835

対数尤度を真数に変換し確率として評価すると

>>> np.exp(model.score(X))
1.3313665376872969e-18

まあ尤度だけ表示したところで「ふーん」って感じなので、次いきましょう。

[復号] 既知のモデルから観測系列 \boldsymbol{x} が得られた時の状態系列 \boldsymbol{s} を復号する

構築したモデルの predict メソッドに観測系列 \boldsymbol{x} を渡すことで、その観測系列が得られる確率が最も高い状態系列が復号される。

上で出力した観測系列サンプルから状態系列を復号してみると

>>> model.predict(X)
array([0, 1, 2, 2, 0, 0, 0, 1, 1, 1])

この結果は上で出力した状態系列と同じなので、尤もらしい状態系列を復号できていることが分かる。

[推定] 未知のモデルから観測結果 \boldsymbol{x} が得られた時のモデルのパラメータを推定する

いよいよパラメータの推定。

新しくパラメータ未設定のモデルを初期化し、fit メソッドに観測系列を与えることでパラメータを推定する。 モデル初期化時の n_iter パラメータによって、推定に使う Baum-Welch アルゴリズムイテレーション回数を指定できる、大きめの値にしておくことで、計算量は増えるが収束しやすくなる。

上で出力した長さ10の観測系列サンプルからパラメータを推定してみると

remodel = hmm.GaussianHMM(n_components=3, covariance_type="full", n_iter=100)
model.fit(X)

結果は以下のとおり (遷移確率と出力値の平均のみ表示)。まあ教師データのサンプル数が10しかないので、まともに推定できていない。

>>> remodel.transmat_  # 推定された遷移確率
array([[  0.00000000e+00,   0.00000000e+00,   1.00000000e+00],
       [  5.00000000e-01,   5.00000000e-01,   5.63659518e-81],
       [  1.99850496e-01,   1.99975083e-01,   6.00174422e-01]])

>>> remodel.means_  # 推定された出力確率 (正規乱数の平均)
array([[  4.43588108,   6.934226  ],
       [ 99.00879995,  99.44400045],
       [  5.87656445,   4.60734161]])

今度は元のモデルからサンプルを10000個出力し、それを用いて推定してみる。

X, Z = model.sample(10000)
remodel = hmm.GaussianHMM(n_components=3, covariance_type="full", n_iter=100)
remodel.fit(X)

今度は、以下のとおり正解に近いパラメータが推定できていることが確認できる。パラメータ内のベクトルの並び順はランダムなので注意。

>>> remodel.transmat_  # 推測された遷移確率
array([[ 0.70077029,  0.1060636 ,  0.19316611],
       [ 0.31953291,  0.39437367,  0.28609342],
       [ 0.2991481 ,  0.19790302,  0.50294888]])

>>> remodel.means_  # 推定された出力確率 (正規乱数の平均)
array([[  2.95450793e-02,   9.44972133e-03],
       [  9.99727677e+01,   1.00005236e+02],
       [  9.99470491e+00,   9.98823474e+00]])

パラメータの初期値を指定して推定する

上の例では初期状態確率 \boldsymbol{\pi}、遷移確率 \boldsymbol{A}、出力確率 \boldsymbol{B} の初期パラメータを指定していないため、すべてのパラメータを乱数で初期化して推定を開始していた。

一方で、パラメータをおおまかに推測できる場合は、パラメータを乱数で初期化する代わりに推定値で初期化することで、収束を早めることができる。

GaussianHMM クラスの init_params パラメータを用いて乱数で初期化するパラメータを頭文字で指定できる。初期値は stmc つまり startprob_, transmat_, means, covars_ のすべてのパラメータを乱数で初期化する設定になっている。よって、例えば大まかな出力の平均値が分かっているというのであれば、次のように init_params パラメータに m を指定せず、means_ プロパティに明示的に初期値を与える。

remodel = hmm.GaussianHMM(n_components=3, covariance_type="full", init_params='stc', n_iter=100)
remodel.means_ = np.array([
        [0.0, 0.0],
        [10.0, 10.0],
        [100.0, 100.0]])
remodel.fit(X)

特定のパラメータを固定し、残りのパラメータを推定する

上の例ではすべてのパラメータを学習プロセスで更新していた。

一方で、一部のパラメータについて正解値が分かっていたり固定値で仮定したい場合には、学習プロセスで更新しないように固定化できる。

GuassianHMM クラスの params パラメータを用いて学習プロセスで更新するパラメータを頭文字で指定できる。初期値は stmc つまりすべてのパラメータを更新する設定になっている。よって、例えば出力の平均値が確定的に分かっているというのであれば、params パラメータおよび init_params パラメータに m を指定せず、かつ means_ プロパティに明示的に初期値を与える。

remodel = hmm.GaussianHMM(n_components=3, covariance_type="full", init_params='stc', params='stc', n_iter=100)
remodel.means_ = np.array([
        [0.0, 0.0],
        [10.0, 10.0],
        [100.0, 100.0]])
remodel.fit(X)