HMMのスケーリング係数

HMMを実装する上で考えなければならない問題点

HMMにおけるパラメータ推定はフォワードバックワードによって効率的に計算される。しかし、実際に計算機上でこのアルゴリズムを動かそうとすると気を付けなければならないことがある。再帰式
\alpha(\mathbf{z}_n) = p(\mathbf{x}_n | \mathbf{z}_n) \sum_{\mathbf{z}_{n-1}} \alpha(\mathbf{z}_{n-1}) p(\mathbf{z}_n | \mathbf{z}_{n-1})
の右辺のそれぞれの確率はよく1より非常に小さい値になり、sequenceの長さが長くなっていくにつれて\alpha(\mathbf{z}_n)が指数的な速度で0に近づいていってしまうという問題がある。こうした場合に計算機上で\alpha(\mathbf{z}_n)をうまく扱えない場合があるので、どうにかしたい*1

i.i.d.なデータ(つまり確率がproductで書ける)なら対数を取ることでこういった問題は避けられたが、この場合は和が入っているのでそう簡単には解決できない。いわゆるlogsumexp*2を使って、これらがきちんと計算できるようにする方法もあるが、ここではscaling係数を使った方法について取り扱う。

scaling係数を使った方法

フォワードバックワードでは\alpha(\mathbf{z}_n)
\alpha(\mathbf{z}_n) = p(\mathbf{x}_1, \cdots, \mathbf{x}_n, \mathbf{z}_n)
と定義した。scaling係数を使った方法では、\alpha(\mathbf{z}_n)の代わりに\hat\alpha(\mathbf{z}_n)を使う。\hat\alpha(\mathbf{z}_n)
\hat\alpha(\mathbf{z}_n) = p(\mathbf{z}_n | \mathbf{x}_1, \cdots, \mathbf{x}_n) = \frac{\alpha(\mathbf{z}_n)}{p(\mathbf{x}_1, \cdots, \mathbf{x}_n)}
と定義される。\alpha(\mathbf{z}_n)ではn + 1個の確率変数の同時分布について考えていたが、\hat\alpha(\mathbf{z}_n)はK個の値を取る1つの確率変数についての確率分布であるから数値計算上の振舞いが\alpha(\mathbf{z}_n)よりよくなるであろう、と推測される。ここで、スケーリング係数c_nを導入する。c_n
c_n = p(\mathbf{x}_n | \mathbf{x}_1, \cdots, \mathbf{x}_{n-1})
で定義される。スケーリング係数を使うと同時分布p(\mathbf{x}_1, \cdots, \mathbf{x}_n)は乗法定理を使ってp(\mathbf{x}_1, \cdots, \mathbf{x}_n) = \prod_{m=1}^n c_mと書き表せる。これを使うと
\alpha(\mathbf{z}_n) = p(\mathbf{x}_1, \cdots, \mathbf{x}_n) p(\mathbf{z}_n | \mathbf{x}_1, \cdots, \mathbf{x}_{n-1}) = \left(\prod_{m=1}^n c_m \right) \hat\alpha(\mathbf{z}_n)
とでき、\alpha(\mathbf{z}_n)\hat\alpha(\mathbf{z}_n)を使った式に書き換えることができた。これを使うと再帰式
\alpha(\mathbf{z}_n) = p(\mathbf{x}_n | \mathbf{z}_n) \sum_{\mathbf{z}_{n-1}} \alpha(\mathbf{z}_{n-1}) p(\mathbf{z}_n | \mathbf{z}_{n-1})
\hat\alpha(\mathbf{z}_n)を使った式に書き換えることができて
c_n \hat\alpha(\mathbf{z}_n) = p(\mathbf{x}_n | \mathbf{z}_n) \sum_{\mathbf{z}_{n-1}} \hat\alpha(\mathbf{z}_{n-1}) p(\mathbf{z}_n | \mathbf{z}_{n-1})
とできる(m個のproductのnのところだけが生き残る)。スケーリング係数c_nをどうやって計算するか?ということに関してはz_{n1}, \cdots, z_{nK}に関して上の右辺を計算して、それらの和の逆数をc_nとすればよい。この係数はあとで使うので記憶しておく。

同様にして、\betaのほうもスケーリング係数を使って書き直すことができて
c_{n+1} \hat\beta(\mathbf{z}_n) = \sum_{\mathbf{z}_{n+1}} \hat\beta(\mathbf{z}_{n+1}) p(\mathbf{x}_{n+1} | \mathbf{z}_{n+1})  p(\mathbf{z}_{n+1} | \mathbf{z}_n)
となる。c_{n+1}\hat\alpha(\mathbf{z}_n)を計算する際に記憶させているので再度計算する必要はない。また、スケーリング係数の積によって、尤度関数を計算することができ、具体的には
p(\mathbf{X}) = \prod_{n=1}^N c_n
とできる(これで前向きアルゴリズム、後ろ向きアルゴリズムで尤度を計算するところはできるようになる...はず)。

また、フォワードバックワードでパラメータ推定するときに必要となる十分統計量のところだが、これも\hat\alpha(\mathbf{z}_n)\hat\beta(\mathbf{z}_n)を使って記述することができる。

各種アルゴリズムのまとめ

コードを書くときに結局どうすればいいんだっけ?となるのでまとめておく。間違っている可能性大。

尤度計算

まず、\hat\alpha(\mathbf{z}_n)を前向きに計算する。

DPを使うので、初期確率\hat\alpha(\mathbf{z}_1)を最初に求める。定義より、\hat\alpha(\mathbf{z}_1) = p(\mathbf{z}_1 | \mathbf{x}_1) = \frac{p(\mathbf{x}_1, \mathbf{z}_1)}{\sum_{\mathbf{z}_1} p(\mathbf{x}_1, \mathbf{z}_1)}とでき、これは初期確率\piとemission parameterを使えば計算できる(\mathbf{x}_1\mathbf{z}_1についての二重ループを作って、zについての和を取るなどすればよい)。c_1 = p(\mathbf{x}_1) = \sum_{\mathbf{z}_1} p(\mathbf{x}_1, \mathbf{z}_1)であるから、\hat\alpha(\mathbf{z}_1)の計算を使ってc_1も計算できる(ので保存しておく)。

alpha_hat = Array.new[Array.new]
c = Array.new

(2..N).each{|n|
  sum = 0.0
  (1..K).each{|k|
    # # (13.59)式の右辺。\mathbf{x}_nはobservedで分かっている
    tmp = p(\mathbf{x}_n | \mathbf{z}_n = k) \sum_{\mathbf{z}_{n-1}} \hat\alpha(\mathbf{z}_{n-1}) p(\mathbf{z}_n = k | \mathbf{z}_{n-1})
    alpha_hat[n - 1][k] = tmp
    sum += tmp
  }
  c[n - 1] = 1.0 / sum
  // 正規化
  (1..K).each{|k|
    alpha_hat[n - 2][k] /= sum
  }
}

scaling係数の積を使って尤度は計算することができる。
p(\mathbf{X}) = \prod_{n=1}^N c_n

パラメータ推定

scaling係数を使った方法でもEMを使うことには変わりない。Mステップでは、Eステップで求めた期待値を使ってtransition parameterなどを更新していくが、ここのステップはscaling係数を使わないときと特に変化はなくてよい。変わるのはEステップ。Eステップで求める期待値が(13.64)式や(13.65)式に変わるのがscaling法で違うところ。

Eステップで期待値を求める前にscaling係数のバージョンの前向きアルゴリズム、後ろ向きアルゴリズムをそれぞれ走らせて、\hat\alpha(\mathbf{z}_n)\hat\beta(\mathbf{z}_n)c_nをそれぞれアップデートしておく。

(あとで書く...)

その他

C++でスケーリング対応のHMMの実装もしています。
www.yasuhisay.info

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

パターン認識と機械学習 下 (ベイズ理論による統計的予測)

  • 作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
  • 出版社/メーカー: 丸善出版
  • 発売日: 2012/02/29
  • メディア: 単行本
  • 購入: 6人 クリック: 14回
  • この商品を含むブログを見る

*1:自分でも実験してみたが、sequenceの長さが大体100を越えた辺りから尤度がinfとかNaNに行ってしまわれた...

*2:logsumexp(log sum exponential)とは - EchizenBlog-Zwei言語処理のための機械学習入門 (自然言語処理シリーズ)を参照のこと