信学会のノンパラメトリックベイズ講座をようやく書き終えました。
具体的な学習例(上のイメージ)や細かい図を描く必要があり, 結局連休後半からずっと かかった気がします。

最後のページではInfinite HMM (NIPS 2001) [pdf] の紹介をしています。 ちょうど岡野原君がohmmをリリースした所で, やたらとタイミングがいいのですが..。
HMMはよく考えるとかなり凄いモデルですが, 上のohmmも含め, 普通のHMMは 隠れ状態の数は事前にセットしておく必要があります。 これに対し, IHMMは隠れ状態の総数すらも観測データを見るだけで決めてくれる という驚異的なモデルで, 僕はD3の時(2003年くらい)に知って, かなり感動しました。

ただ, IHMMは理論を理解するのもそうですが, 実装がかなりややこしいので *1 僕は実際に実装はしていなかったのですが, 最近素晴らしいことに, Ghahramaniの ところにいる Van Gael がIHMMのMATLABツールキットを 公開 してくれているので今回それを使ってみました。 以下, ドキュメントに書かれていないことを含めて, 自分用も兼ねて使い方のメモ。

• ihmmでは, 観測値 y は「行ベクトル」で与えます。つまり, y = [23 41 57 8 ..] のような整数の列ということ。これは y = load('y.txt')'; のようにして作れます。
• 実際の推定は

  >> [states,stats] = iHmmSampleGibbs(y,hypers,200,1,1,ceil(rand(1,length(y))*2));
  Iteration 0: K = 2, alpha0 = 0.317225, gamma = 0.742441.
  Iteration: 1: K = 2, alpha0 = 0.129917, gamma = 1.106346, JL = -68126.516533.
  Iteration: 2: K = 2, alpha0 = 0.234067, gamma = 0.560162, JL = -67684.190698.
  Iteration: 3: K = 2, alpha0 = 0.108580, gamma = 1.092227, JL = -67309.148779.
  Iteration: 4: K = 2, alpha0 = 0.333919, gamma = 0.634966, JL = -67124.438335.
  Iteration: 5: K = 2, alpha0 = 0.201366, gamma = 0.676564, JL = -66761.077048.
  :
  

  のように行います。上は2個の隠れクラスから学習を始める例。 この時, ハイパーパラメータの入った構造体 hypers をyと一緒に 渡す必要がありますが, これには H, alpha0_a, alpha0_b, gamma_a, gamma_b という フィールドが必要で, 特にHには, emission確率 p(y|k) のディリクレ事前分布の ハイパーパラメータを渡す必要があります。
  その場でいちいち指定すると面倒なので,

  function hypers = mkparam(y)
  % hypers = mkparam(y)
  % returns default hyperparameters in ihmm.
  L = max(y);
  hypers.H = ones(1,L) * 0.2;
  hypers.alpha0_a = 4;
  hypers.alpha0_b = 2;
  hypers.gamma_a = 3;
  hypers.gamma_b = 6;
  

  のような関数を作って, hypers = mkparam(y); とすると便利です。
• 重要: Gaelのツールキットは, MATLABなのであまり大きなデータを想定していないため, yがあまり大きいと(5000〜10000以上), 途中で謎の所で落ちます。
  これは何故かというと, 内部で使っている第1種 Stirling 数のテーブルを超えて しまうからなので, 基本的には付属の Stirling1.m を使って

  s = Stirling1(2000,2000);
  save(s, 'StirlingCache.mat');
  

  としてキャッシュを増やしておけばOKです.. *2 が, スターリング数は組み合わせ的な数で, 簡単に言うと竹内関数みたいなもの なので(そこまでではありませんが), Stirling1(10000,10000); などとすると 計算が終わらなくて死ねます。
• そこで最初, どちらにしてもスターリング数の数が大きい場合は近似が効くはず なので, 近似すればいいんじゃね?と思って少し調べたところ, Hwang氏の "Asymptotic Expansions for the Stirling Numbers of the First Kind" [citeseer] などの文献が見つかり, なるほどと思っていました。
• ところが, 実際にスターリング数のマトリックスを見てみると, 当然ですが大きい 数の場合は, 行列の内部は境界を除いてほとんど Inf です。 たとえば, s(2000,2000) = 1 ですが,
  • s(2000,1999) = -1999000
  • s(2000,1990) = 2.7247e+56
  • s(2000,1950) = 1.5842e+250
  • s(2000,1900) = Inf
  です。よって, まともにスターリング数を計算する必要はほとんどなく, 実際には 次のような関数 snumbers1.m を用意して,

  function a = snumbers1(n,m)
  persistent snumbers
  if isempty(snumbers)
    snumbers = load('StirlingCache.mat');
    snumbers = abs(snumbers.s);
  end
  K = size(snumbers,1);
  
  if (n == m)
    a = 1;
  else 
    if (n > K) | (m > K)
      a = Inf;
    else
      a = snumbers(n,m);
    end
  end
  

  iHmmHyperSample.m の snumbers(N(j,k),l) を snumbers1(N(j,k),l) に置き換えて やると, 無事大きなデータでも動くようになります。

というわけで, 上のようにすれば, 観測値yを用意するだけで, Infinite HMMを 実際に試すことができます。上では離散値(例えば単語のID)ですが, iHmmNormal* という関数もあるので, 正規分布でもOKのようです。

ちなみに, van Gaelはもちろん Beam Sampler を用意してくれていて iHmmSampleBeam()を使うとビームサンプリング (状態数が無限の場合でも可能な, Slice samplingを使ったForward-Backward) ができますが, 僕が試した限りでは普通のGibbsの方が計算や収束が速かったようです。 系列がこのツールでは一つだけなので, Gibbsの方が隠れ状態を増減する機会が多く, 柔軟なのかもしれません。

---