mots quotidiens. | |
Daichi Mochihashi (持橋大地) daichi <at> ism.ac.jp | by hns, version 2.10-pl1. |
|
||||||||||||||||||||||||||||||||||||||||||||||||||
θ = [0.4, 0.3, 0.2, 0.1] のような離散分布をランダムに初期化したいと いうことは, 自然言語処理や混合モデルの学習でよくある状況だと思う。 下で書くようにこれはガンマ分布からのサンプリングに還元できるので, MCMCなどのベイズ学習一般にもよくある問題。
さて, θは適当に [0,1] の一様乱数で初期化してもいいのだが, 値がかなりバラバラに なってしまうので, 例えば [0.2609, 0.2836, 0.1974, 0.2581] のように 「ある値を中心としてそこから少しずれた」ように初期化したい時は, θ ~ Dir(α) とディリクレ分布からサンプリングすればよい。 ディリクレ分布 Dir([α1,α2,..,αK])からのサンプルを取るには, ガンマ分布に従う独立なサンプルγk ~ Ga(αk, 1) (k = 1 .. K)を発生させて, それを和が1になるように正規化して
θk = γk / Σiγiとすればよい。ここでΓ分布の確率密度関数は
Ga(x|a,b) = a^b / Γ(a) x^(a-1) exp(-bx).MATLABには gamrnd() という関数があるので, MATLABでディリクレ分布からの サンプルを生成するには
function theta = dirmult(alpha) % theta = dirmult(alpha) % Multinomial sampler from a Dirichlet distribution. % alpha : row vector of Dirichlet parameters % theta : row vector of Multinomial output % $Id: dirmult.m,v 1.1 2004/09/14 04:39:40 dmochiha Exp $ theta = normalize(repmap(alpha,@gamrnd_plus,1)); function x = gamrnd_plus(alpha,beta) if alpha == 0 x = 0; else x = gamrnd(alpha,beta); endのようなコードを書けば, 簡単にできる。(α_i = 0 の場合を特別扱いしているのに注意。)
これは1年以上前に書いたコードなのだが, Cにはガンマ分布からの乱数を生成
するような便利な関数はないので, どうしようか, ということが問題になる。
簡単な場合には
Nealが紹介している方法
があるが, a が実数の場合や, 1よりずっと小さい場合(自然言語ではこれが普通)
では使えない。
MATLABでは gamrnd はどう実装されているかと思って調べてみると
(Statistics Toolbox が入っている場合は, 'type gamrnd' とタイプしてみるべし),
Devroye (1986) "Non-Uniform Random Variate Generation", Springer-Verlag
に書かれている方法を使っているらしい。
あれ, これはちょうど大羽さんの
ベイズWiki
で最近伊庭さんが
紹介された
濃い本では..(笑)ということで,
このページ
から本が全部PDFでダウンロードできるようです。
これを見て, Cで実装してみたものが下です。
/* random.c $Id: random.c,v 1.1 2005/12/27 04:11:12 dmochiha Exp $ References: [1] L. Devroye, "Non-Uniform Random Variate Generation", Springer-Verlag, 1986 http://cgm.cs.mcgill.ca/~luc/rnbookindex.html */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include "random.h" double gamrand (double a) { double x, y, z; double u, v, w, b, c, e; int accept = 0; if (a < 1) { /* Johnk's generator. Devroye (1986) p.418 */ e = exprand(); do { x = pow(RANDOM, 1 / a); y = pow(RANDOM, 1 / (1 - a)); } while (x + y > 1); return (e * x / (x + y)); } else { /* Best's rejection algorithm. Devroye (1986) p.410 */ b = a - 1; c = 3 * a - 0.75; do { /* generate */ u = RANDOM; v = RANDOM; w = u * (1 - u); y = sqrt(c / w) * (u - 0.5); x = b + y; if (x >= 0) { z = 64 * w * w * w * v * v; if (z <= 1 - (2 * y * y) / x) { accept = 1; } else { if (log(z) < 2 * (b * log(x / b) - y)) accept = 1; } } } while (accept != 1); return x; } } void dirrand (double *theta, double *alpha, int k, double prec) { int i; double z = 0; /* theta must have been allocated */ for (i = 0; i < k; i++) if (prec != 0) theta[i] = gamrand(alpha[i] * prec); else theta[i] = gamrand(alpha[i]); for (i = 0; i < k; i++) z += theta[i]; for (i = 0; i < k; i++) theta[i] /= z; } double exprand (void) { return (- log(RANDOM)); }"充分に発達した科学は、魔法と区別がつかない" という言葉がありますが, なんかほとんど黒魔術みたいなコードですね。(笑)
#define RANDOM ((double)rand()/(double)RAND_MAX)とするのが簡単です。(これは工藤君のコードからコピーしたような気がします。 Thanks.) MCMCを動かすなどで正確な乱数が必要な場合には, MTなどに 置き換えればよいかと思います。 これは基本的にMATLABで実装されているものと同じで, 実験してみましたが 平均, 分散やヒストグラムもほぼ同じでした。とは言っても, 一応無保証という ことで。
ディリクレ分布からサンプルする場合は, 上の dirrand を使って, double *theta, *alpha をアロケートしてから,
dirrand(theta, alpha, k, 0);
とすれば theta に Dir(alpha) からのサンプルが得られます。 最後の引数はディリクレ分布を中心 α/Σ(α) と精度 Σ(α) に分解してサンプリング する時に, 精度(中心への集中の度合い)を調整できるようにするためのパラメータ.
タイトル一覧 |