2012-03-16

スライスサンプリングでお手軽サンプリング

ベイズ法などで用いるサンプリング法のなかで，スライスサンプリングというのがあります．原論文は以下になります．
SLICE SAMPLING (Neal, Annals of Statistics 2003)
サンプリング法で代表的なのはMetropolis-HastingやGibbs Samplingですが，これらにはそれぞれ問題点があります．例えばMetropolis法の問題点として，実際の問題に適応する場合，どのような提案分布を選んだら良いのか分からない，ということが挙げられます．よく使われるのは，現在のサンプル点を中心としたGauss分布ですが，分散をどう選べば良いか，ということなど，自明ではありません．提案分布が適切でない場合，サンプルの棄却率が高くなり，アルゴリズムの効率が落ちます．一方，この分散を小さくしすぎると，一度のサンプルで動く距離が小さくなり，これも効率の低下につながります．またGibbs Samplingの場合，サンプルしたい一変数を，他の変数で条件付けた分布からのサンプルを考えるわけですが，この分布が簡単に求まらない場合もあります．
このような問題に対処するため，スライスサンプリングは，提案分布を自分で考えなくても，現在のサンプル点に応じて自動でサンプル幅（上の分散に対応）を選んでくれるサンプリング法です．つまり，難しいことを考えなくてもどんな問題にも適応可能なわけです．そういう意味で，お手軽なサンプリング法と言えるでしょう．
さて，このスライスサンプリングですが，今回簡単に試せるライブラリを見つけたので紹介します．C++のコードになりますが，Mark JohnsonのAdaptor Grammarのパッケージの中に入っています．
Software available from Mark Johnson, Dept of Computing, Macquarie University, Sydney, Australia

実験

ちょっと簡単に試してみました．
メインで呼び出すメソッドは，

template <typename F, typename LogF, typename Uniform01>
F slice_sampler1d(const LogF& logF0,               //!< log of function to sample
                  F x,                             //!< starting point
                  Uniform01& u01,                  //!< uniform [0,1) random number generator
                  F min_x = -std::numeric_limits<F>::infinity(),  //!< minimum value of support
                  F max_x = std::numeric_limits<F>::infinity(),   //!< maximum value of support
                  F w = 0.0,                       //!< guess at initial width
                  unsigned nsamples=1,             //!< number of samples to draw
                  unsigned max_nfeval=200)         //!< max number of function evaluations

です．logF0というのは，求めたい分布の，正規化されていないxの関数値の，対数を返す関数オブジェクトです．xには，1つ前のステップで求めた（つまり現在の）サンプル点を，u01には，[0,1)を返す乱数生成オブジェクト，を入れます．求めたい分布の関数で，対数を指定しないといけないのは，まあそういう決まりということで・・・．

例えば，適当な比率になっている2つの混合Gauss分布からのサンプルを求めるコードは，以下のようになります．

#include <iostream>
#include <iomanip>
#include "slice-sampler.h"

using namespace std;

struct LogGaussianMixture {
  LogGaussianMixture(double mean1, double mean2, double variance, double pi)
      : mean1(mean1), mean2(mean2), variance(variance), pi(pi) {}
  double operator()(double x) const {
    double y1 = exp(-(x - mean1) * (x - mean1) / (2 * variance)) + 1e-100;
    double y2 = exp(-(x - mean2) * (x - mean2) / (2 * variance)) + 1e-100;
    return log(pi * y1 + (1 - pi) * y2);
  }
  double mean1;
  double mean2;
  double variance;
  double pi;  
};
struct RandGen {
  double operator()() const {
    return (double)rand() / RAND_MAX;
  }
};

int main(int argc, char *argv[])
{
  LogGaussianMixture log_gm(0, 10, 3, 0.2);
  RandGen rand_gen;

  double x = 0;
  int sample_size = 1000000;
  for (int i = 0; i < sample_size; ++i) {
    x = slice_sampler1d(log_gm, x, rand_gen);
    cout << setprecision(5) << x << " ";
  }
  
  return 0;
}

ちょっと注意なのが，関数値にlogをかける前に，微小な値を足しています．これは，関数値が0になった場合に，log(0) = -∞となり，ライブラリ中のassertでこけるためですｗ．
出力したファイルをRで描画してみたのが，以下になります．

具体的なアルゴリズム

これまでサンプリングの中身には全く触れませんでしたが，あまりにブラックボックスだと信用出来ないと思うので，簡単に解説します．
以下，一次元の場合で考えますが，多次元でも同じです．
今，ある分布p(x)がサンプルをとりたい分布だとして，p(x)が正規化されていない，関数f(x)∝p(x)が手元にあるとします．これはベイズで事後分布を求めるときに自然に出てくる話です．
このとき，スライスサンプリングは，以下の手順でp(x)からサンプリングを行います．

x を適当に初期化
以下を繰り返す
1. 現在のxの値x0での，関数の値f(x0)を計算する．
2. [0, f(x0))の一様分布から，値yをサンプリングする．
3. そのyの値に対して，y <= f(x) となる，xの範囲Sを求める．（これは，現在のyの値で関数f(x)を水平に"スライス"し，f(x)と交わる点を求めることに相当）
4. 範囲Sの一様分布からxをサンプルし，これを新しいxとする．

これは結局何をしているかというと，y = f(x)をグラフに書いて，f(x)より下にある部分の点(x,y)を，移動させるわけです．この(x,y)のうち，xだけを取り出すと，確率f(x)に沿った値でサンプルが得られる，という話になります．
ここで問題になるのが，y <= f(x)となるxの範囲Sを求める部分で，これはf(x)の関数形が複雑な場合自明ではありません．これに対処するため，論文中には2つのアプローチが紹介されています．両方の共通点は，最初に現在のサンプル点であるx0の回りの小さい区間を考え，それを何らかの枠組みで，Sを近似するほど大きく（大きすぎた場合は小さく）していくということです．Mark Johnsonの実装では，このうち，区間を（確率的に左右どちらかの方向に）倍に倍に増やしていく，doubling procedureというのに基づいています．

上の議論は，もともとの論文中での主張で，主に連続分布を対象にしています．一方，離散分布の場合でもスライスサンプリングは別の意味で有効性が示されていますが，こちらは[機械学習] スライスサンプリング - tsubosakaの日記などを参照してください．

もちろん実際の問題に対しては，もっと良いサンプリング法が個別に存在することもあると思いますが，とりあえず手軽に試せるということで，サンプルしたいパラメータの分布が（正規化されてなくても）書き下せるときは，応急処置的にこれを使って対処してみる，というのは良いのではないかと思います．このパッケージは一次元のサンプリングしか対応していませんが，多次元の場合，自明な方法はGibbs Samplingに組み込む（GibbsもNLPなどでよく出てくる自然な分布だと簡単に定式化出来ますが，一変数のサンプルが自明でない場合もある）ことで対応できます．また，ハイパーパラメータのサンプルなどであれば，多くの場合一次元ですが，尤度関数が非常に複雑ということが多いため，このサンプラーが有効だと思います．

2011-08-12

A Hierarchical Pitman-Yor HMM for Unsupervised Part of Speech Induction

論文紹介です。今日引っ越す関係で、しばらく家からネットが繋がらなくなりそうで、自分の理解の整理も含めて書いてみます。まとまってないです。
今年のACLで、A Hierarchical Pitman-Yor Process HMM for Unsupervised Part of Speech Induction [ pdf ]というのがある。これは、教師なしの品詞推定モデルで、やったのはBlunsomという、最近だとTree Substituion GrammarをPYPでモデル化する？ということとかやった人。ジャーナルは読もう読もうと思って読めてない状態なので適当ですが。この論文は、教師なし品詞推定に、過去の系列モデル（言語モデルや品詞推定）のモデル化におけるstate-of-the-artな手法を様々組み合わせたら、過去最高の結果が出たよ、というもの。この論文はACLが公開されてすぐに読んだけれど、その時は推定手法を飛ばしてざっと読んで、ふーんという感じだった。しかしもう一度読んでみて、結構推定のところで工夫しているのがこの論文の肝みたいなので、それについてだらだらと書いてみる。
このモデルの新しいところは、（多分）HMMの階層化を考えたところ。つまり、３グラムHMMの推定に２グラムとの補間を用いる。すごいシンプルな気がするけど、これが自然に扱えるようになったのがPYPやDPで階層化が自然に出来るようになったから、だと思う。
で、これの推定で何が問題になるかというと、階層CRPにおけるテーブルのカウントが問題になる。推定は他のBayesian HMMと同じくGibbsを用いるのだが、このとき３グラムをサンプルしている途中で、サンプルの候補によってその後の確率が変化するということが生じる。（この辺りはちゃんと説明するとかなり面倒なので割愛します。といってもこれを説明しないと後の議論は成り立たないのですが・・・。余裕ができたら補足するかも）で、本来ならば３グラムをサンプリングしている途中で、テーブルが増えるかどうかをサンプリングしながら推定するのだろうが、それは重すぎるので、計算途中では近似的にやってしまおう、ということをしていて、それが式（３）になる。
まず最初に言っておくと、この式の分母は多分タイポで、E_n[K_i]b^Uとなっているところのb^Uは、a^Uの間違いのはず。そう考えると、これは現在のテーブル数の近似値に、新しくテーブルが増える確率を足して、新しい近似値とする、という形になっていて、ある程度納得がいく。で、これが実際のPYPのテーブルの期待値とどれぐらい差が出るのかというのが示されている。式（３）のこともあるので、この図を簡単にシミュレートしてみた。
C++のコードが以下のような感じ。

#include <iostream>
#include <vector>
#define RANDOM ((double)rand()/(double)RAND_MAX)

using namespace std;
vector<double> approxTableCount(int customer, double a, double b) {
  double p_0 = 0.25;
  vector<double> approxCounts(customer, 0);
  
  for (int i = 0; i < customer; ++i) {
    double approxE = i == 0 ? 0 : approxCounts[i-1];
    approxCounts[i] = approxE + (a * approxE + b) * p_0 / ((i - approxE * a) + (a * approxE + b) * p_0);
  }
  return approxCounts;
}

vector<double> expectTableCount(int customer, double a, double b) {
  double p_0 = 0.25;
  vector<double> expectCounts(customer, 0);
  for (int j = 0; j < 1000; ++j) {
    int tableCount = 0;
    for (int i = 0; i < customer; ++i) {
      double random = RANDOM;
      if (random < (a * tableCount + b) * p_0 / ((i - tableCount * a) + (a * tableCount + b) * p_0)) {
        ++tableCount;
      }
      expectCounts[i] += (tableCount - expectCounts[i]) / (j+1);
    }
  }
  return expectCounts;
}

int main(int argc, char *argv[])
{
  vector<double> approxTables = approxTableCount(100, 0.9, 1);
  vector<double> expectTables = expectTableCount(100, 0.9, 1);
  for (size_t i = 0; i < approxTables.size(); ++i) {
    cout << (i+1) << " " << approxTables[i] << " " << expectTables[i] << endl;
  }
  return 0;
}

この結果をプロットしてみると、以下のようになる

このように、図3の a = 0.9の場合が再現出来ているので、式（３）は間違いだったことになる。
論文ではこのあと、1 tag per word の制限、つまり、１つの単語は複数の品詞を持たない、という制限を置いてこの近似を使って推定すると、精度が上がるといっている。具体的には、サンプルする単位が、各typeに対応する品詞となり、関係する同時確率を計算し、比較してサンプリングするのだと思う。最初はこれがtype based sampling（Liang 2010）と関係あるのかと思って、訳が分からない、と思っていたのだが、どうももっとシンプルな話で、単にその語に関係するトライグラムを全て抜き出して、その語の品詞が k であるとした時の全同時確率を計算し、サンプリングする、という話になるよう（ちゃんと理解出来てない気がするので、間違ってたら教えて下さい）。

あと教師なし品詞推定について思うことを少し。これも含めて、ここ１〜２年の教師なし品詞は、1 tag per word制限のもとで求めるというのが主流で、そっちのほうが精度が良い、ということらしい。この精度が良いというのは、既存の品詞体系と、ここで求まった品詞体系とを比較して、最もよくマッチするものが取れているということで言っている。しかし、品詞というのはそれを求めることが目的にはならないもの（自然言語処理的には、係り受け解析や翻訳の前段階）で、教師なしでgold standardに近づけて意味があるのかは微妙な気がする。まして、1 tag per wordの制約を入れてしまうと、精度は上がるかもしれないが、本質的な品詞とは何かという問題には何も答えないんじゃないかと思う。この辺りは、murawakiさんが似たことを書かれている。
だけど、少し視点を変えると、これは完全にテキストデータのみから品詞っぽいものを取り出す方法を示していて、これは教師あり学習が行えるほどのリソースが存在しない言語の解析に役に立つんじゃないかという気がしている。というのも、最近インターンで、言語横断的な何かをするということをしていて、リソースの少ない言語にどうアプローチするかというのが、今後問題になりそうな感じなので。リソースのない言語に対して、どうにか情報を取り出したいという場合、精度はそこそこでもそれっぽいものが取れていれば良くて、そういう問題には1 tag per wordも有効に働くのかなあとも思う。教師なしのモデル化というのは、個人的には人間の思考モデル、言語の本質を捉えるための道具（うまくはまれば、人間をそういうモデルで近似出来る）、と考えていたのだけれど、工学的意義を考えると、これはこれで意味があるのかという視点で見ることが出来て、面白いなあと思ったり。

2011-07-18

distance dependent Chinese Restaurant Process

お久しぶりです。nokunoさんにも紹介されてしまったので頑張って月1ぐらいは更新したいと思ってます…。今回は面白かった論文の紹介です。去年のICMLのBleiの論文で、相変わらずCRPとかです。ICML版はこちらで、longer versionもあり、こちらからダウンロード出来ます。ICML版でほとんどの部分は説明されてて、理論的に詳しいところが知りたい人はlonger versionも補足的に読むといいかもしれません。
以前のエントリーで、DPMを説明するときに、CRPを介して説明出来るということを書きました。これはつまり、データをクラスタリングする場合、データの事前分布にCRPを仮定し、CRPの事後分布（レストランの状態の分布）がどうなるかを考え、同じテーブルに座ったデータ=客を同じクラスタとすることでクラスタリング出来ることを意味しています。この場合、背後にあるDPという構造を考えずに、単にCRPを道具としてクラスタ数を制限しないクラスタリングを考えています。distance dependent CRP（以下dd-CRP）はこの考えを推し進め、もはやDPとの等価性は存在しないが、より柔軟なクラスタリングを可能にしようというモデルです。
CRPでは、一人の客=データは、どこかのテーブル=クラスタに対応付けられました。dd-CRPでは、一人の客は、もう一人の別の客に対応付けられます。つまり、CRPではテーブル1に客{1,2,3}が座っているという状況は、dd-CRPでは例えば、テーブル1 <- 客1 <- 客2 <- 客3、というリンク構造で表されます。ここで矢印で繋がっている客が1つのクラスタを形成します。なんのためにこんなことをするかというと、例えば系列データをクラスタリングする場合、距離が近いデータは同じクラスタに含まれやすいだろうという直感を入れることが出来ます。例えばニュース記事などが時系列に並んでいてクラスタリングする場合、同じ時期の記事は同じクラスタになりやすい、という予想を組み込むことが出来ます。これはdd-CRPの生成過程において、新しい客が行く先の客の分布に対して、近くの客ほど行きやすい、というような分布を仮定したり、行く先の客との距離に閾値を設けたりすることで達成されます（論文ではdecay function）。ここで面白いのが、この閾値を定めず、今より前の客のところであればどこでも同じ確率で行く、とすると、これは従来のCRPと同じものになります。ただdd-CRP自体は、見て分かるように、客=データの位置関係によって同時確率が定まるため、これは交換可能な列とはなりません。つまり背後にDPは仮定出来ないということで、ノンパラベイズの枠組みとも異なる、と書かれています。
この事後分布はGibbsで推定出来ます。データを抜き取り、客の再配置を行うわけですが、このときその客には、他の客からのリンクが貼られている可能性があります。つまり5 <- {3,6,8 <- {2}}という具合。この場合{2,3,5,6,8}は同じクラスタを形成し、5に3,6,8がくっつき、8が2にくっついています。ここで5の再配置を行う場合、それが従えるデータの集合{2,3,5,6,8}を全て一緒に再配置する、ということになります。これが従来のCRPだと、5の再配置を行う場合、5が座るテーブルを変えるだけで、他の{2,3,6,8}については変えません。そのため、dd-CRPを従来のCRPの置き換えとして使ってこの推定を行うと、より速く混合が進み、マルコフ連鎖の収束が早まるらしいです。他の実験では（これがメインですが）、一応新聞記事や時系列に並んだNIPSデータセットで、このモデルが従来のCRPよりも低いパープレキシティを示すと書いてありました。
個人的にはこの収束が早まるというのがこのモデルの面白いところかなと思いました。現状dd-CRPは単なるCRP、DPMのモデルの置き換えとして使えるようですが、現在主に使われるのは階層モデル（HDP）だと思うので、HDPについてもdd-CRPのような等価な表現が出てくれば、色々嬉しいことがありそうです。特に隠れ変数が強い相関を持つモデルの場合従来のCRP、Gibbsでは収束が遅いことがボトルネックなので、この辺が改善されれば面白そうです。

ちなみに、Blei本人がこれのプログラム（R）を公開しています。
http://www.cs.princeton.edu/~blei/downloads/ddcrp.tgz
ちょっとだけ触ってみましたが、従来のCRPとの収束の違いなどはシミュレート出来ないようです。（いくつかパッケージが入っていることが前提です）また、僕の場合safelogという関数が定義されていない、と怒られたのですが、よく分からないのでここにあるsafeLogという関数を取ってきて名前をsafelogに変えました（いいのか？）。

> library(plyr)
> library(lda)
> library(Matrix)

# coraはldaに入っている論文のabstractを集めたデータセット
> data(cora.document) # 各論文毎にBOW形式で入っている
> data(cora.vocab)    # スカラと単語との対応データ

> source('data-modeling.R')
> source('ddcrp-inference.R')

> dat <- corpus.to.matrix(cora.documents[1:100],cora.vocab) # 文章データをBOWの疎行列に変換
> res <- ddcrp.gibbs(dat=dat[1:100,], dist.fn=seq.dist, alpha=10, # run gibbs
                     decay.fn=exp.decay(5),
                     doc.lhood.fn(0.5), 10, summary.fn = ncomp.summary)

> res$summary # クラスタ数の遷移を表示
        iter ncomp
summary    0   100
           1    10
           2     8
           3     9
           4    11
           5     9
           6     7
           7     7
           8     5
           9     7
          10     6
          11     6
          12     6
          13     7
          14     6
          15     6
          16     6
          17     7
          18     8
          19     9
          20     6
          
> res$map.state # 最後の状態を表示
    idx cluster customer
1     1       1        1
2     2       1        1
3     3       3        3
4     4       1        1
5     5       1        2
6     6       1        1
7     7       1        4
8     8       1        2
9     9       1        7
10   10       1        9
11   11       1        7
12   12       1       10
13   13       1        6
14   14       1       11
15   15       1       10
...
47   47       1       42
48   48      16       45
49   49       1       47
50   50      50       50
51   51       1       47
52   52       1       51
53   53       1       51
54   54       1       53
55   55       1       52
56   56      16       46
57   57       1       54
58   58       1       57
59   59      16       56
60   60       1       49
61   61      16       59
62   62      16       61
63   63       1       54
64   64      16       61
65   65       1       63
66   66       1       60
67   67      16       61
68   68       1       66
69   69       1       63
70   70       1       58
...

まあ、これだけじゃよく分からないですが…。最後の出力は、最後のレストランの状態を示していて、左からデータ=客のインデックス、その客が属するクラスタ番号、対応する（くっついている）客のインデックス、を示しています。idxとcustomerを見ると、近くの客のところに行きやすいようになっていることが分かります。
ちなみにcoraのデータセットを使ってるのは、使えるのがこれしか見当たらなかったからで、これは時系列データではないのでその点で微妙です。推定のddcrp.gibbs()では、いくつか引数を指定します。decay.fnは、上で述べたdecay functionで、ここではexp.decayという指数的に遠くの客に行きにくくなるものを使っています。decay.fn=window.decay(dim(docs)[1])とすると、前に存在する全ての客に等しく遷移するという、従来のCRPと等価なdd-CRPのサンプルが得られます。
Rは久々に触りましたが、パッケージが色々なソフトに分散している状況は何とかならないのですかね。まあ現状matlabもRもnumpyもほとんど触っていない（使えない）のですが。自分で書くのはどれか一つでいいにしろ、どの言語もある程度使えるようになっておく必要性というのは感じます。

2011-06-25

iHMMと単語分割と

何となく復習がてら、Tehによる階層ノンパラベイズのチュートリアル pdf を眺めてたのですが、ちょっと気になったことがあったので、まとめてみます。
infinite Hidden Markov Model（iHMM）というのがあります。名前から推測出来そうですが、これは通常の隠れマルコフモデルの状態数をデータに決定させてしまおうというモデルで、本質的な部分は前回紹介したDPMと全く同じです。これはどういう仕組みかというと、HMMの状態を無限（データから推測）とするために、HMMの状態の集合を、DPからのサンプルとします。DPは加算無限次元の離散分布を出力するので、これでHMMの状態をデータに適応させるような状態が作り出せることになります。具体的には、以下のようなグラフィカルモデルで表現出来ます。（上の論文から借用）

ここで、
$\90dpi G_0|\gamma, H \sim DP(\gamma, H)$
$\90dpi G_\theta | \alpha, G_0 \sim DP(\alpha, G_0)$
となります。 $\90dpi \theta_k$ が隠れ変数で、そこから $\90dpi x_k$ が出力される、というモデルです。で、この $\90dpi \theta_k$ は前の状態 $\90dpi \theta_{k-1}$ によって定まる $\90dpi G_{\theta_{k-1}}$ から出力されます。つまりこの $\90dpi G_{\theta_{k-1}}$ が、通常のHMMにおける確率テーブル $\90dpi P(\theta_k|\theta_{k-1})$ に相当します。で、このテーブルの要素がデータに適応し自動で増えていくように、無限個まで考えるわけですが、HMMにおいて、遷移候補の状態の集合は、どの状態でも共有しなければなりません。例えば状態3からは1,2,3,4にいけるが、4からは1,2,3にしかいけないということになると、HMMでなくなります。これをきちんと定めるために、各 $\90dpi \theta$ で定まる $\90dpi G_\theta$ の基底分布には、共通の離散分布を置かないといけません。DPにおいて基底分布を離散分布とおくと、そのサンプルは、基底分布が持つ値と同じ値のみに確率を持つ、別の離散分布が得られます。例えば{a:0.5, b:0.3, c:0.2}という離散分布があったとします。これを基底分布としてDPのサンプルをとると、例えば{a:0.8, b:0,1, c:0.1}などが得られるわけです。そして、この離散分布の要素数を、加算無限個にしたいのですが、これをするために、 $\90dpi G_\theta$ の基底分布 $\90dpi G_0$ が更に別のDPから生成された、とします。この基底分布 $\90dpi H$ は連続のものをおきます。こうすることで、加算無限個の要素をHMMの状態間で共有出来るモデルを構成することが出来ます。階層モデルにおける推定の詳細などは、(Teh ACL2006) pdf が（比較的）平易です。iHMMの推定はヘビーでビームサンプリングが有効らしいのですが、(Gael EMNLP2009) pdf は、それを用いて教師なしで品詞数をデータから決めてクラスタリングするというものです。

生成モデルによる単語分割

以上が前置きでして、上の論文のiHMMの応用のところで、単語分割の問題が挙がっています。教師なしの生成モデルによる単語分割は、日本では持橋さんの教師なし形態素解析 pdf が有名ですが、もとはSharon Goldwaterの音素列からの単語の獲得、という問題がベースになっています（はず）。詳細は彼女のD論 pdf にあるのですが、赤ちゃんが言葉を習得するときに、連続音声から単語区切りをどう学習しているか、をモデル化する問題です。英語だと単語がスペースで区切られているので、この問題は音声の区切りでしか出てこないのですが、日本語だと書き言葉も連続してるので、以下は日本語の単語分割に置き換えて考えます。
この問題はざっとですが春頃勉強していたのですが、考えてみるとiHMMの一種と言えるなあというのが、今回の主題です。その視点は今までなかったので。この問題では、単語分割されていない文章から、単語の列を得るわけですが、隠れ状態を各単語と考えると、単語の種類は実質無限個ですので、iHMMで考えることが出来ます。持橋さんのモデルは単語トライグラムを階層Pitman-Yor過程でスムージングするというモデルで、詳細は省きますが、要はその時点で（文章を分割して）得られた単語から、次の単語への遷移（次にどこで区切るか）を考えることになります。この遷移確率は、各単語を文字の並びから考えた生成モデルを基底分布とする、PYPから得られます。PYPはPitman-Yor Processの略で、DPの特別バージョンです。これは上で述べたiHMMにおいて、 $\90dpi G_0$ を別のDPからのサンプルで考えていることと実質同じです。違うのは、ここでは $\90dpi H$ を単純な離散分布とするのではなく、文字の組み合わせや長さなどを考慮した単語の生成モデルを考えることで、より精密にしていることです。また $\90dpi \theta$ が各単語に相当し、 $\90dpi G_\theta$ が次の単語の生成確率を定めます。持橋さんのモデルでは単語トライグラムを考えるので、実際には前とその前の語により状態が定まるのですが、これは通常のHMMの次元数を大きくすることと同じです。
ちなみにノンパラのHMM関連ということで、今年のACLで階層PYP-HMMというモデルを考え、品詞を教師なしで推定するもの pdf がありますが、これは基底分布にstaticな離散分布を置くシンプルなモデルです。つまり品詞数を事前に決めてしまうというもので、iHMMとは異なります。

持橋さんはこれをNested PYPと呼んでます。PYPの基底分布に無限の要素が出力可能な別のPYPを置くことで、無限個の単語の出力を可能にした単語Nグラムの生成モデルを考えています。Nested〜は一般的な名称らしく、他の論文ではNCRPが使われています。正確な定義は勉強不足なのですが、NDP、NPYPなどは、iHMMの特別な場合と見なせるのかもしれません。持橋さんの生成モデルが理解しにくかった人は、単語が隠れ状態であるHMMであるとしてみると、違った理解が得られるかもしれません。

2011-06-06

DPMから学ぶノンパラベイズの思想

はじめまして。そろそろ何かしら情報を発信していく必要性を感じたため、主に研究関連で、まとまったことがあれば記事にしていくことにしました。どれだけ更新出来るかは謎ですが。今回は、ノンパラベイズの基本をディリクレ過程を中心にまとめます。

機械学習におけるノンパラベイズは、出て来てから10年以上経っていることもあり、大分一般的な話題になってる気がしますが、例えばブログできちんと分かりやすく説明したものってほとんどないように思います。僕がそもそも研究系のブログをあまりチェックしないというのもあるかもしれないですが、、、。個人的には去年の夏頃からの卒論で、Tehや持橋さんなどの論文を泣きながら読みつつ理解出来なかったので、その時の気持ちを思い出しながら書いてみたいと思います。例えばディリクレ過程（以下DP）を理解しようとして論文など読むと、DPはCRPと等価であるとか、SBPと等価であるとか書いてあって、そもそもDPって何？という状態になるんですよね（僕だけ？）。そんな人の疑問を払拭出来たらなと。個人的な理解も含まれるので、間違っている箇所があったら申し訳ないです。ちなみに、全くディリクレ過程を知らないという人は他の人の記事や持橋さんのスライド[1]などを見るといいかもしれません。またベイズについてPRML1章程度の知識は前提とします。ノンパラとは、（パラメタの個数=モデルの複雑さ）に制限をかけないという意味です。

ノンパラメトリック ベイズ

ノンパラベイズの応用範囲は多岐に渡るため、ここでは最も基本的かつ（思想的に）重要なディリクレ過程混合モデル（DPM）に焦点を当てます。DPMを用いる主な動機は、混合モデルによるクラスタリングにおいて、コンポーネント数kをデータから自動的に決めたい、というものです。例えば2次元データを混合正規分布に従ってクラスタリングしようとした場合、従来kはモデルのパラメタとして与えられるとしていました。kが問題になる場合、モデル選択（周辺尤度の最大化や、交差検証）に基づきモデルの複雑さ=kを決定します。この問題を、kが観測データに応じて自動的に調整されるようなモデルを置くことで、つまりkの推定自体をモデルに組み込むことで自然に解決しよう、というのが、ノンパラベイズの基本的な発想になります。

モデルの設計

さてこのような前提のもと、どうしてDPが出現するのか、DPとは何なのか、について説明します。上で述べたモデルの生成過程は、次のように書けます。
$\90dpi x_i \sim F(\theta_i)$
$\90dpi \theta_i \sim G$
ここで、 $\90dpi F(\theta_i)$ は各コンポーネントの分布で、例えば正規分布。 $\90dpi \theta_i$ がそのパラメタに相当します。我々が推定したいのは、 $\90dpi \theta_i$ が従う分布Gです。これは各 $\90dpi \theta_i$ を要素として持つような離散分布です。従来のパラメトリックな手法では、このGを、各データがどのクラスタから発生するかを規定するインディケータに関する多項分布と、そのクラスタの実際のパラメタに関する分布によって表現します。ベイズの枠組みでは、インディケータである多項分布に対する事前分布としてディリクレ分布を置いて推定したりします。この場合、ディリクレ分布の事後分布を推定出来れば、多項分布のパラメタが分布の形で得られます。*1。ノンパラベイズではこの考えをより一般化します。つまり、クラスタ数kを決めた上で何らかの事前分布（例えばディリクレ分布）を設計するのではなく、kがいくつになるか、まで含めた確率モデルというものを考えます。具体的には以下のCRPでkの個数まで含めた分布を考えることが出来ます。

Chinese Restaurant Process (CRP)

CRPは以下のような生成過程です。

レストランに無限のテーブルが存在し、客が一人ずつ入店する
新しく来た客は、他の客が多くいるテーブルに座りやすい（客が多い = そのテーブルの料理はおいしいのだろうと推測する）
客は、誰も座っていない新しいテーブルに座ることも出来る。この場合、そのテーブルに新しい料理が運ばれる。１つのテーブルにはそれぞれ１つの料理のみが対応する。

CRPは割と総称的な言葉で、具体的なモデルによってバリエーションがありますが、共通部分を抜き出すと上のようになります。よくある図ですが、イメージは下のような感じです。これはレストランに５人の客が座っていて、次の客がどこに座るかを表す図です。それぞれのテーブルに、そのテーブルに座っている客の数に比例する確率で座りますが、CRPのパラメタである $\90dpi \alpha$ に比例する確率で、新しいテーブルに座ります。

これをデータの生成過程と考えます。つまり、新たな客がテーブルに座る毎に、それに対応する料理がデータとして生成され、観測されるというモデルを考えます。このとき重要なのは、CRPが交換可能（exchangeable）な過程であるという点です。交換可能というのは、系列データ $\90dpi x_1, x_2, \cdots, x_n$ に関して、その順序を入れ替えた列 $\90dpi x_{\sigma(1)}, x_{\sigma(2)}, \cdots, x_{\sigma(n)}$ を考えると、この２つの列の生成確率が等しいことをいいます。例えば $\90dpi x_1, \cdots, x_n = 1 2 1 3$ という列が得られたとして、この順序を変えた列 $\90dpi x_1, \cdots, x_n = 3 2 1 1$ というものを考えても、その同時確率は変わりません。これはCRPを実際にシミュレートして計算してみるとすぐに分かると思います。実際、同時確率を式で表せば証明出来ます。
ここで、レストランの客を実際に観測される各データとすると、同じテーブルに座っている客を、同じクラスタに属するデータと考えることが出来ます。つまり１つのテーブルが１つのクラスタに相当し、そこに座っている客が、そのクラスタに属するデータ、とするわけです。そして、各テーブルに対応する料理は、そのクラスタに対応する正規分布のパラメタ（平均分散）に対応付けます。このように考えると、上で述べたような、データに応じてクラスタ数が自動的に増加するようなモデルの生成過程を定義出来たことになります。ここで問題なのは、CRPは逐次的な生成過程であるということです。つまり、１つ１つのデータがそれまで生成されたデータ（に対応する客がどのテーブルに座ったか）に応じて、生成されるというモデルです。我々が推定したいのは、混合モデルの混合比、及び各コンポーネントのパラメタを定める１つの離散分布が何かしら存在し、その分布に応じてコンポーネントが選択され、各データが生成される、というモデルです。イメージは以下のような感じです。

実はこのモデルをBayes推定した結果とCRPとが一致するという結論になるのですが、これを理解するには、数学的な定理としてde Finettiの表現定理というものを用います。これは以下のようなものです。
ある系列 $\90dpi x_1,\cdots,x_n$ が交換可能な性質を持つとする。このとき、この結果はあるパラメタのもとでの独立同分布な試行の結果として表現することが出来る。そのパラメタを積分消去した確率と、交換可能な生成過程で生成されたとする同時確率は等しくなる。
$\90dpi p(x_1,\cdots,x_n) = \int p(\theta) \prod_{i=1}^{n} p(x_i | \theta) d\theta$
結論をいうと、この $\90dpi p(\theta)$ がディリクレ過程と呼ばれるものです。これは数学的な結果なのですが、上の定理において、パラメタである $\90dpi \theta$ は分布（測度）に置き換えることが可能です。これを明示的に $\90dpi \theta$ の代わりにGを用いて書くと
$\90dpi p(x_1,\cdots, x_n) = \int p(G) \prod_{i=1}^{n} p(x_i | G) dG$
となります。上で説明したようにCRPの生成列は交換可能なので、この $\90dpi x_1,\cdots ,x_n$ をCRPで生成されるデータとおけます。さて右辺ですが、まずこのGの１つは具体的な離散分布を定めます。データはGから独立同分布に得られ、それをp(G)のもとで平均をとるという式です。これはまさに、離散分布Gに対する事前分布p(G)を考え、そのもとでGを周辺化（積分消去）するというベイズの予測の式です。ベイズ的には、データを発生させた真の分布（パラメトリックな問題では分布のパラメタ）を１つ決めるのでなく、その分布が従う分布を求めます。その上で、新たなデータに関する予測としては、データを発生させたとする分布を、積分消去したものが欲しいわけです（つまりp(G)自体は必ずしも必要ではない）。このp(G)が、ディリクレ過程と呼ばれるものです。そしてp(G)からのサンプルが、いわゆる加算無限次元の離散分布といわれるものです。要するに、我々は本当は上の式の右辺を計算したかったが、p(G)をどのように設計すればいいのかよく分からない。代わりにCRPが交換可能であることを用いて、CRPから左辺で計算する、ということが出来るのです。さらに推定の際には、Gibbs samplingで交換可能であることが本質的に重要になります。この辺の詳細は他の論文などに譲ります。代表的なのは[2]など。

周辺の話題

またDPの説明の際Polya Urn表現というのも出てきますが、これはCRPと全く同じです。Stick Breaking Process（SBP）についてもよく出てきますが、上では説明しませんでした。これはSBPがなくてもDPについて説明出来るからです。一応解説すると、SBPはDPから具体的な１つのサンプル分布Gを得る方法です。つまり、SBPから具体的なGが構成出来ますが、このサンプルを取りまくれば、p(G)が原理的に得られるというわけです。推定の際は、CRPを基にすると、交換可能性からGibbs samplingを用い、SBPを基にすると変分ベイズを用いて推定出来ます。

まとめ

要するに解釈の仕方なのですが、歴史的な経緯はどうであれ、DPが機械学習のどういう要請で出現したのかが重要なのだと思います。繰り返しになりますが、DPは混合モデルにおける混合数を実質無限まで拡大出来るようにしたい、という要請から生まれたものです。でも実際どんな分布を置けばそれを表現出来るかというと難しい、、、。代わりにde Finettiの定理を足がかりに、無限個まで要素を増やせる生成モデルを考えてみればいいのでは？となります。それが交換可能な系列を生成するのであれば、de Finettiの定理から、そのような分布の存在を暗に証明することが出来るのです。
以上の解釈でノンパラベイズを捉えると、理解しやすいと思うのですがどうでしょうか。例えば、Dirichlet Diffusion Tree（DDT）[3]というのがあるのですが、これはDPMを一般化したモデルです。DPMの各クラスタが、それぞれ独立ではなく階層的な構造を持っているとするモデルで、ノンパラモデルの一種です。これに関しても、データの背後に、階層モデルを生成する１つの真の木構造があると仮定します。そしてベイズ推定で、その木を積分消去したものを考えるのですが、これは、データを階層的に生成するCRPのような交換可能な生成過程を考えることで、上の議論と全く同様に、真の木構造を暗に仮定したモデルを考えることが出来ます。DPの場合は一応SBPによって具体的なサンプルを直接得ることが可能ですが、このような複雑なモデルになると直接具体的なモデルのサンプルを構成することは困難です（存在を暗に仮定するだけ）。DDTに関しては今年のICMLでGhahramaniらが推定を高速化したらしく、今個人的にチェックしてます[4]。

参考資料

[1]Nonparametric Bayes for Non-Bayesians, 持橋大地 [IBIS2008] pdf
[2]Markov Chain Sampling Methods for Dirichlet Process Mixture Models, Neal, R.M. [2000] pdf
[3]Density Modeling and Clustering Using Dirichlet Diffusion Trees, Neal, R.M. [2003] pdf
[4]Message Passing Algorithms for Dirichlet Diffusion Trees, Knowles, D.A. et al. [ICML 2011] pdf

*1:この辺りは専門外なのですが、このように混合正規分布をベイズ推定する話はあまり見ないので一般的ではないのかもしれません。言語処理的には、オリジナルのLDAが混合モデルの混合数（トピック数）を決めた上でベイズ推定する代表でしょうか。しっくりこない人は適当に読み替えながら進んでください。