mots quotidiens. | |
Daichi Mochihashi (持橋大地) daichi <at> ism.ac.jp | by hns, version 2.10-pl1. |
|
|||||||||||||||||||||||||||||||||||||||||||||||||
いま, 離散的なデータ
θ = {θ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) に等しい,ということにバスの中で気付いた。
Γ(α+n) / Γ(α).ここで α は実数で, n は整数(単語のカウント)。 実際には上の値は n が大きくなるとかなり大きな値になるので, 対数を取った lnΓ(α+n) - lnΓ(α) が必要になる。
よって, 原理的にはこの比を計算するためには Γ関数は必要なく,
log(α) + log(α+1) + .. + log(α+n-1) を計算すればよい。
ただし, n が大きくなるとこの方がかえって計算が重い可能性があるので,
α(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 が小さいほど掛け算が少ないので, 速い。
もちろん, この高速化は, そもそも 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(x) は(glibcでは)どうやって実装されているのだろうか, と思って調べてみたら, Intelのプロセッサには fyl2xp1 という命令 (解説) があって, それを使っているらしい。何と。つーか知らないよ, そんなこと。(笑)
ディレクターに色々質問をして話を聞いた後,
工藤君と形態素解析の話をしたり, 高林君と夕食を食べたりして帰りました。
ちなみに,
去年のこの時期
には, 広尾の統数研に行っていたようです。
タイトル一覧 |