ぐはっ。第1種Stirling数..。
どうしてここでそれが出てくるんだ。ぐはっ。(吐血)

やっと理解できた。何日も悩んでいたので, 自分用メモというか, 公開デバッグ。

いま, 離散的なデータ θ = {θ₁, ..., θ_n}が Dirichlet Process G ~ DP(α_0,G_0) に従っていたとする。
つまり, θ ~ G, i.i.d. n times.
このとき, n個のデータθが, 実際には k 個の違った値をとる確率 p(Z_n = k) が知りたい。(Z_n は n 個のデータのもつカテゴリの個数を表す 確率変数.)
例えば, θ = {{θ₁,θ₂},{θ₃,θ₄,θ₅},{θ₆,θ₇}}とグループ分けされるとき, その確率は (α_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(Z_n = k) = Σ_i∈π(n,k) α^k・f(i)/α⁽ⁿ⁾ = α^k/α⁽ⁿ⁾Σ_i∈π(n,k)f(i)= α^k/α⁽ⁿ⁾ 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=1ⁿ s(n,k) α^k/α⁽ⁿ⁾.
これを上の式と比べると,
p(Z_n = k) = s(n,k) α^k/α⁽ⁿ⁾
となることがわかる。□
直接考えても, (Σ_i∈π(n,k)f(i))は (n個のデータをk個に分割する組み合わせの数) × (各partitionの中での円順列の数) なので, 第1種Stirling数の定義に従って これは s(n,k) に等しい,ということにバスの中で気付いた。

論文ばかり読んでいて煮つまってきたので。
DM等の計算では, 次のようなΓ関数の比を計算する必要がよく生じる。

Γ(α+n) / Γ(α).

ここで α は実数で, n は整数(単語のカウント)。 実際には上の値は n が大きくなるとかなり大きな値になるので, 対数を取った lnΓ(α+n) - lnΓ(α) が必要になる。
前はこの式をナイーブに gammaln(α + n) - gammaln(α) として計算していたのだが, Γ関数には Γ(x+1) = xΓ(x) という性質が あるので(部分積分すると出てきたはず),
Γ(x+n)/Γ(x) = (x+n-1)Γ(x+n-1)/Γ(x) = ... = (x+n-1)‥xΓ(x)/Γ(x) = x(x+1) .. (x+n-1) となって, 下の日記の記号を使うと, これは raising factorial x⁽ⁿ⁾ に等しい。

よって, 原理的にはこの比を計算するためには Γ関数は必要なく, log(α) + log(α+1) + .. + log(α+n-1) を計算すればよい。
ただし, n が大きくなるとこの方がかえって計算が重い可能性があるので, α⁽ⁿ⁾が int に収まる範囲(αが小さい値の場合, n < 12くらい)は log(α(α+1)..(α+n-1)) を返して, n がそれ以上大きい場合には gammaln(α+n) - gammaln(α) を返すように すれば計算が速くなるはず。
というわけで, 計算してみた。

% ./g 2.5 1 1000000
gammaln(2.5 + 1) - gammaln(2.5), iteration = 1000000 times
gammaln:    0.513084 seconds.
gammaratln: 0.181330 seconds.
% ./g 2.5 12 1000000
gammaln(2.5 + 12) - gammaln(2.5), iteration = 1000000 times
gammaln:    0.630065 seconds.
gammaratln: 0.313247 seconds.

n が小さいほど掛け算が少ないので, 速い。
ということで, Γ関数を明示的に計算しないことで, だいたい2〜3倍は速くなる ようだ。

もちろん, この高速化は, そもそも lnΓ(x) = gammaln(x) がどのくらい複雑な計算に なるかに依存している。 gammaln(x) は一番簡単には, 下のような感じで書けるようだ。 (gamma(x) は, exp(gammaln(x)) とすればよい。)

/*
 *  gamma.c
 *  Lanczos expansion of the gamma function.
 *  from "Numerical Recipes", Section 6.1.
 *  $Id$
 *
 */
#include <math.h>
#include "gamma.h"

double
gammaln (double x)
{
        static double coeff[6] = {
                76.18009172947146,
               -86.50532032941677,
                24.01409824083091,
                -1.231739572450155,
                 0.1208650973866179e-2,
                -0.5395239384953e-5
        };
        double y, z, s = 1.000000000190015;
        int i;

        y = x;
        z = x + 5.5;
        z -= (x + 0.5) * log(z);
        for (i = 0; i < 6; i++)
                s += coeff[i] / ++y;
        return (- z + log(2.5066282746310005 * s / x));

}

これは "Numerical Recipes in C" のプログラムを冗長な所を除いて, わかりやすく書き直したもの。 実際には上の実験には Lightspeed に含まれている実装を使っていてもう少し複雑だが, これでも, x < 400 くらいの場合は誤差は 10^(-5) 以下 くらいにはなるようだ。この簡単な版を使うと, 速度差は最大で4倍くらいになる。
ただ, x^(n) がオーバーフローしないようにする等の考慮が大変なことを 考えると, 数倍程度の高速化だと割に合わないかもしれない。 intensive に上の比を計算する場合は役に立ちますよ, ということかも。

＠ log:

高林君の招待で, 東京の Google ラボに行って話してきました。
30分でプレゼンして下さい, とのことで話したのだけれど, はっきり言って かなり難しかったと思うので惨敗。orz
「確率分布の確率分布」という概念(ディリクレ分布のことね)を出さずに 話そうと思ったのですが, それは無理だったことが判明。;
ちなみに Dirichlet Process は,「無限次元の多項分布の確率分布」という ものなので, これもまた専門でない人には視覚的に説明しにくい概念のような 気が..。
専門でない人向けに説明するには, もっと別の直感的なスライドを作った方が よかったかも。(と言っても, 言語モデルの導入や, 前に伊藤君にちょっと話した 応用等を書くだけで6時間とかその位かかっていて, 30分のために全部書き直す 余裕はなかったのだけれども。)
それでも, 出た質問はかなり鋭いものが多く, 終わった後にさすがだな と思いました。

ディレクターに色々質問をして話を聞いた後, 工藤君と形態素解析の話をしたり, 高林君と夕食を食べたりして帰りました。
ちなみに, 去年のこの時期 には, 広尾の統数研に行っていたようです。