今日のぞくパタ読書会の予習に「続・わかりやすいパターン認識」(以降「ぞくパタ」)の 11章をつらつら読む。
p227 に「CRP における使用テーブル数の変化」というグラフがあり、αが 2 と 10 のそれぞれにおいて、来客数が 1000 になるまでシミュレーションしたときのテーブル数の推移が示されているのだが、当然ながら一回の試行しか描かれていない。
中華料理店過程(CRP)ないしホップの壺 をちょっこり試すだけなら 10行ちょっとのスクリプトでできる話であり、実際やってみると、実行するたびにテーブルの増え方が結構変わる。つまり一回の試行ではぶっちゃけよくわからん。
というわけで、α=2 について 1000 回シミュレーションを行い、来客数がそれぞれ 100, 200, 400, 600, 800, 1000 になったときのテーブル数の分布を箱ひげで描いてみた。
ほとんどは中心の周辺にきれいにまとまるけど、散らばりもそれなりにある。
また、テーブル数がある程度増えた以降は、分散はほとんど増えなくなっていく（散らばり具合は n が少ないうちにだいたい決まっちゃう）ということもわかる。テーブルの増える確率がどんどん小さくなることから想像できる通りではあるが。
ちなみに、それぞれの平均値と標準偏差はこんな感じだった。
100: mean=8.31(8.39) sd=2.47(2.42) 200: mean=9.69(9.77) sd=2.78(2.68) 400: mean=11.03(11.14) sd=3.07(2.93) 600: mean=11.86(11.95) sd=3.20(3.06) 800: mean=12.41(12.53) sd=3.25(3.15) 1000: mean=12.88(12.97) sd=3.27(3.22)
カッコ内はそれぞれの理論値。1000 サンプルにしては悪くない？
平均値(期待値)は前回の記事「ディリクレ過程(中華料理店過程)のトピック数(テーブル数)の期待値を導出してみる #ぞくパタ - Mi manca qualche giovedi`?」で求めていたが、分散も「独立な確率変数の和の分散は、分散の和」*1であることに注意すればまったく同様に求めることができる。
ここでは厳密さが欲しかったわけじゃないのでシミュレーションでお茶を濁したけど、実は二手間くらいかければ確率質量関数も計算できる。もし知りたかったら「ポアソン二項分布」で検索。
大したものではないが、シミュレーションのスクリプトも一応貼っておく。
中華料理店過程は、事後分布とか絡まなければ、実装はとても易しく楽しいので、興味あればこんなん見ないで自分で実装してみる方がいい。
import numpy alpha = 2.0 labels = [100, 200, 400, 600, 800, 1000] target = dict() data = [] for i, x in enumerate(labels): target[x] = i data.append([]) for i in xrange(1000): n = [] for N in xrange(1001): z = N + alpha p = [ni / z for ni in n] p.append(alpha / z) k = numpy.random.choice(len(p), p=p) if k<len(n): n[k] += 1 else: n.append(1) if N in target: data[target[N]].append(len(n)) for i, N in enumerate(labels): p = alpha / (numpy.arange(N) + alpha) m = p.sum() v = ((1 - p) * p).sum() print "%d: mean=%.2f(%.2f) sd=%.2f(%.2f)" % (N, numpy.mean(data[i]), m, numpy.std(data[i]), numpy.sqrt(v)) import matplotlib.pyplot as plt plt.boxplot(data, labels=labels) plt.show()
*1:正確には無相関で十分
