こちらもどうぞ。
numpy は R と同じように出来ると書いたけど、特にループを出来るだけ廃したければ、いろいろコツが必要。
しばらく放っておいたら忘れそうなので、メモ。
R で数式を実装する場合にも似たような考察は必要なので、参考になるかも？
隠れマルコフ with EM アルゴリズムの数式を１つずつピックアップして、それを実装するには、というチュートリアルっぽく。
まずは隠れマルコフモデルで真っ先に実装する式がこれ(Baum-Welch のフォワードステップの初項)。
一見複雑な数式だが、実は下の簡単なベクトルを求めればよいことがわかる。わか……らないと、実装できない(苦笑)*1。
z が 1-of-K、つまり (0,..,0,1,0,..,0) という K 個の要素の内どれか１つだけ１であるベクトルであることに注意。
これを計算するには、πと p(x|φ_k) が必要。
πは多項分布(1-of-K)を与えるものなので、足して1になる K 個の値。これはベクトルで与えるべき。初期値は均等もしくはランダム。
均等なら 1/K が K 個並んだベクトルを作る。numpy には同じ値を繰り返してベクトルを作る、というメソッドが見あたらないので、 全部1のベクトルを作る numpy.ones() を使おう。
numpy.array([1/K] * K) でもいいが、ones() なら複数次元でも全く同じように使えるのでちょっと嬉しい*2。
import numpy K = 6 # なにか適切な値が入っている pi = numpy.ones(K) / K # 1/K が K 個並んだベクトル
【追記】
numpy.tile() というのがあった。複数次元もOK。numpy.repeat() ってのもあり。
pi = numpy.tile(1.0 / K, K) # 1/K が K 個並んだベクトル
numpy.ones() の方がタイプ数少ないけどね。
【/追記】
ランダムに初期化するなら、ディリクレ分布からサンプリングすれば自動的に「足して1になる」ので便利。
pi = numpy.random.dirichlet([1] * K) # 足して1になる K 次ベクトル(ランダム)
ディリクレ分布に与えるパラメータは欲しい乱数の傾向を反映する。
端っこの値はあんまり嬉しくない(1/K にノイズが載るくらいでいい)なら (10,10,...,10) とかもっと大きい値にすればいいし、わからない(無情報)なら (1,1,...,1) でいい。
ちなみに EMA-HMM の話だが、πの初期値を Dir(0.1, ..., 0.1) からサンプリングすると、pi のいくつかの要素がつぶれるが、尤度はおおむね上がる。まだどういう理屈でそうなるのかよくわかっていない(事前分布じゃあなくて単なる初期値なのになあ)。
p(x|φ_k) は HMM の出力確率であり、モデルによって全然違うわけだが、今は離散分布(given トピックに対する単語の生起確率)を考える。
単語の種類数が V 個とすると、１つのトピックに対し多項分布(1-of-V)になるので、πと同じ要領で V 次元のベクトルになる。これがトピックの数だけあるので、V×K の2次元の量になる。
R だと迷わず行列で持てばいいが、Python + numpy の場合は numpy.array(ndarray) で持つか、numpy.matrix で持つか、Python の list(ネイティブな配列)を使うか、という選択肢が出てきてしまって、迷う。
行列として積や逆行列を取る場合は matrix、内包表記( R の apply や sapply みたいなの)を使いたい場合は list、列ベクトルを取りたいケースがあれば ndarray を使う、というのが目安かな。ここでは ndarray でいこう。
このモデルでの p(x|φ_k) は対称な値で初期化しておけば良さそうだから、全ての値が 1/V である V×K の2次元 arrayを作ろう。
V = 300 # 何か正しい値が入っている # emission transition B = numpy.ones((V, K)) / V # 全ての値が 1/V である V×K の2次元 array
ndarray の積 * は要素ごとの掛け算なので、 は2つのベクトルの積を取るだけでよい。
alpha = [] # alpha の系列を格納する list x = [1,5,6,9] # 入力系列(単語IDの列) N = len(x) # 入力系列の長さ alpha.append( pi * B[x[0]] ) # alpha(z_1) = (.., pi_k p(x_1|φ_k), ...)
alpha は N×K の量になるのだが、K 次のベクトルを list に格納している。実装が一番短そうだと思ったので。あとで内包表記使うし。
もし ndarray に入れたいなら、
alpha = numpy.zeros((N, K)) alpha[0] = pi * B[x[0]]
とすればいい。
そうそう、R と Python + numpy のもう一つ大きな違いがあった。
配列／ベクトル／行列のインデックスが R では 1 から始まるが、Python + numpy は 0 から始まる。
数式(よっぽどのことでもない限り 1 から始まる)をコードに落とすとき、R はそのままでいいが、Python の場合は読み替えないといけないので、うっかり間違えないようにしないといけない(そして間違ってしまったら、そのバグは見つけにくい……)。
って、たった1行の数式の説明でこんなにかかるのか。
本当は Baum-Welch のξの更新式とか、Viterbi の max とるところとかがめんどくさかったから、そこらへんをメモしたかったのになあ。企画倒れ？
というわけで、続きは気が向いたら。
