やっと理解できた。何日も悩んでいたので, 自分用メモというか, 公開デバッグ。
いま, 離散的なデータ
θ = {θ1, ..., θn}が Dirichlet Process
G ~ DP(α_0,G_0) に従っていたとする。
つまり, θ ~ G, i.i.d. n times.
このとき, n個のデータθが, 実際には k 個の違った値をとる確率
p(Z_n = k) が知りたい。(Z_n は n 個のデータのもつカテゴリの個数を表す
確率変数.)
例えば,
θ = {{θ1,θ2},{θ3,θ4,θ5},{θ6,θ7}}とグループ分けされるとき, その確率は (α_0 = αと書くと)
1・1/(α+1)・α/(α+2)・1/(α+3)・2/(α+4)・α/(α+5)・1/(α+6)
= 2・α^3/(α(α+1)‥(α+6)) = α^3・1^(2)/α^(7)
になる。ここで, x^(n) は raising factorial x(x+1)‥(x+n-1).
同様にして, n 個のデータが k 個に分かれるとすると, その partition のパターン
を π(n,k) と書くと,
p(Zn = k) = Σi∈π(n,k) αk・f(i)/α(n) = αk/α(n)Σi∈π(n,k)f(i)= αk/α(n) c(n,k)
となる。ここで c(n,k) は n と k から決まる定数。
いっぽう, 第1種Stirling数 s(n,k)
(reference)
(Mathworld)
を使うと, raising factorial x^(n) は
x^(n) = s(n,1) x + s(n,2) x^2 + .. + s(n,n) x^n
と展開できるので, 両辺を x^(n) で割って x = α とおくと,
1 = Σk=1n s(n,k) αk/α(n).
これを上の式と比べると,
p(Z_n = k) = s(n,k) αk/α(n)
となることがわかる。□
直接考えても,
(Σi∈π(n,k)f(i))は (n個のデータをk個に分割する組み合わせの数) ×
(各partitionの中での円順列の数) なので, 第1種Stirling数の定義に従って
これは s(n,k) に等しい,ということにバスの中で気付いた。