3月の 言語処理学会大会2006 の チュートリアル で、筑波大学の山本先生と二人で講師をすることになりました。
基本的に山本先生がされるので, 僕は最後の30分か15分かで話すことになるかと 思います。
そうすると, 講師のメンバーの中で松本研は坪井君, 工藤君, 僕, 乾君 と4人もいることになります。 そもそもチュートリアル資料の提出先が三人とも松本研関係者だったりして, 松本研の寡占状態というか, 下で書いた Power law distribution だなあと思ってしまいます。
(もちろん, あくまで国内に限った話ですが。)

本当は「意味に基づく統計的言語モデル」とかいうタイトルにしたいところ ですが, 言語処理学会の言語学方面で言う論理式的なものとは少し違うので, こんなタイトルになりました。

Bayesian NLP workshop については, 主催者の Hal Daume III が blog で summary を書いているようです。(下で書いた SIGGRAPH の話も含まれています。) 僕も少し思ったコメントを書いていたりして。
別のエントリ (NLP and NIPS) には, NLP関連の今年のNIPSの論文(ベイジアンに限らない)が まとめられているようです。

やった。見事に最初の方のコンポーネントしか使われていないのに注意。 しかし, 大きなデータで試すとダメだったりする。; バグがー。

θ = [0.4, 0.3, 0.2, 0.1] のような離散分布をランダムに初期化したいと いうことは, 自然言語処理や混合モデルの学習でよくある状況だと思う。 下で書くようにこれはガンマ分布からのサンプリングに還元できるので, MCMCなどのベイズ学習一般にもよくある問題。

さて, θは適当に [0,1] の一様乱数で初期化してもいいのだが, 値がかなりバラバラに なってしまうので, 例えば [0.2609, 0.2836, 0.1974, 0.2581] のように 「ある値を中心としてそこから少しずれた」ように初期化したい時は, θ ~ Dir(α) とディリクレ分布からサンプリングすればよい。 ディリクレ分布 Dir([α₁,α₂,..,α_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));
}

"充分に発達した科学は、魔法と区別がつかない" という言葉がありますが, なんかほとんど黒魔術みたいなコードですね。(笑)
Γ分布のような複雑な分布の場合は, rejection を使ってサンプルを生成する ようです。 ここで RANDOM は [0,1] の間の一様乱数を返すマクロで,

#define RANDOM ((double)rand()/(double)RAND_MAX)

とするのが簡単です。(これは工藤君のコードからコピーしたような気がします。 Thanks.) MCMCを動かすなどで正確な乱数が必要な場合には, MTなどに 置き換えればよいかと思います。 これは基本的にMATLABで実装されているものと同じで, 実験してみましたが 平均, 分散やヒストグラムもほぼ同じでした。とは言っても, 一応無保証という ことで。
お持ち帰り用パッケージ: random.h random.c

ディリクレ分布からサンプルする場合は, 上の dirrand を使って, double *theta, *alpha をアロケートしてから,

dirrand(theta, alpha, k, 0);

とすれば theta に Dir(alpha) からのサンプルが得られます。 最後の引数はディリクレ分布を中心 α/Σ(α) と精度 Σ(α) に分解してサンプリング する時に, 精度(中心への集中の度合い)を調整できるようにするためのパラメータ.

-