論文ばかり読んでいて煮つまってきたので。
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分のために全部書き直す 余裕はなかったのだけれども。)
それでも, 出た質問はかなり鋭いものが多く, 終わった後にさすがだな と思いました。

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