Chinese Restaurant Process を発生させてみた

  • ????????????????????

タイトルの通りです。単なる遊びです。Blackwell-McQueen Polya’s urn scheme とか表記されたりもするみたいです。

Chinese Restaurant Process (CRP)

よく(名前にもあるように)中国料理屋の例が出ますが、個人的には「初期の壺に連続的なスペクトルの色(虹色)を持つボールが入っているバージョンのポリアの壺」と考えるほうがわかりやすかったです。つまり、以前に記事にしたポリアの壺の色数が(可算とは限らない)無限集合に拡張されたもの、ということです。兎に角、CRP では以下のスキームで無限の長さの確率変数列を発生させます。

  1. \(G_0\) を適当な確率分布、\(\alpha\) を正の実数とします
  2. \(k\leftarrow1\)とする
  3. \(G_0\) から乱数(色)をひとつサンプリング
  4. その乱数を記録し、同じ色のボールを空の壺にいれる
  5. \(k\leftarrow k+1\)とする
  6. \(\alpha/(\alpha+k-1)\) の確率で \(G_0\) からサンプリングし、それ以外のときは壺からボールを一つ取り出してそれをサンプルとする
  7. サンプルと同じ色のボールを壺に入れる
  8. ステップ5に戻る

実際には無限に数列をつくることはできないので、10000個サンプルした段階でやめることにします。長さ10000の系列を15回生成して、その分布をプロットしてみます。\(G_0\) は標準正規分布、\(\alpha=3\) としました。

# Generating n random variables via CRP(alpha, G0)
rcrp<-function(n,alpha,G0){
  res<-numeric(n); res[1]<-G0(1);
  for(i in 2:n){
    p<-alpha/(alpha+i-1)
    res[i]<-ifelse(runif(1)<p,G0(1),sample(res[1:(i-1)],1))
  }
  class(res)<-"crp"
  res
}

# Simple simulation
par(mfrow=c(5,3)); par(mar=rep(1.5,4)); par(oma=c(1,1,1,1))

N <- 10000

for(i in 1:15){
  plot(table(rcrp(N,3,rnorm))/N,xlim=c(-3,3),ylim=c(0,0.4),ylab="",xaxt="n")
  axis(side=1,at=-2:2)
  curve(dnorm,add=T)
}

ここからわかること:

  • \(G_0\) は連続値をとる密度関数だが、得られる分布は離散的(実はこれは10000個で生成をやめずに、無限個生成しても離散的)
  • 10000個のサンプルのほとんどは数~十数のポイントに落ちる

いろいろ調べると(証明すると)わかること:

  • 実は得られた15個の離散分布はディリクレ過程(←ここでは説明していない)に従う「ランダム確率分布」である
  • (これは、もちろん有限で止めているという意味では正確な表現ではない)
  • 本当のディリクレ過程からのサンプルされる分布は(離散だけども)サポートは無限集合になる
  • CRPによって得られる確率変数列は過去に依存しているので、iidには見えない。しかし、実際にはこれらは得られた離散分布からのiid とみなすことができる
はてなブックマーク - Chinese Restaurant Process を発生させてみた
Pocket