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$ をそれぞれアップデートしておく。

(あとで書く...)