stats(統計学会メーリングリスト)やibismlに流れましたが, 来月1/19(水)の 統数研NOE統計的機械学習セミナーで話すことになりました。 [Link]

前半は Yee Whye で後半が僕という感じで, 僕は今やっている半教師あり形態素解析 の話をする予定です。
この話は外で発表するのは今回が最初です。 今年の言語処理学会2011 にも出す予定ですが(きのう申し込んだ), NLP2011は豊橋技科大なので, 東京方面でご興味のある方はどうぞお越し下さい。

Unsupervised and Semi-supervised Learning of Nonparametric Bayesian
Word Segmentation

For the second part, I extend the NPYLM to semi-supervised learning using
Conditional Random Fields (CRF). Although NPYLM can be regarded as a kind of
semi-Markov model, naive combination with semi-Markov CRF is prohibitive
and proves to work badly. To cope with this problem, we convert the
information between Markov CRF and semi-Markov NPYLM to yield a consistent
combination of discriminative and generative models.  We show the results
on segmenting twitters, speech transcripts, dialects based solely on
newspaper supervised data, as well as the results for standard datasets on
Chinese word segmentation.

* Latter half of the talk is a joint work with Jun Suzuki and Akinori Fujino
(NTT CS Labs).

学習プログラムを動かす際に, パラメータ設定その他が後で わからなくなってしまうことがあるので(注: 自然言語処理はデータ量が多いので, 学習は通常数時間から数10時間かかる), それら+αの記録のために, logging というクラスを定義しておくことを 思いついた。 *1

class logging {
        logging (char *argv[], char *file);
        ~logging ();
};

のようなクラスを定義しておくと, mainで

#include "log.h"
int main (int argc, char *argv[])
{
	while ((c = getopt(argc, argv, opts)) != -1)
		何とか..
	logging start (argv, model);
	学習本体;
}

のようにするだけで, model.log にコンストラクタでargv,起動時刻,ホスト名 *2 などの情報が 書き込まれ, デストラクタで終了時刻が自動的に記録される。以下のような感じ。

init   : ./spylm -n 1 -e 1e-8 -a 4 -b 1 -N 2 ../data/genji.txt genji.model
host   : samurai5.cslab.kecl.ntt.co.jp
start  : Mon Dec 29 00:06:56 2008
--------
iter[0] char height = 9
iter[0] char height = 9
--------
finish : Mon Dec 29 00:10:39 2008

start というのはただの名前ですが, こうやって使うといかにも記録をスタート したかのように見えると思います。
logging.cpp の中身は大体こんな感じ。

#include <cstdio>
#include <cstdlib>
#include <cstdarg>
#include <unistd.h>
#include <ctime>

static FILE *lp = NULL;

logging::logging (char *argv[], char *file)
{
        char *s, hostname[BUFSIZ];
        time_t timer;

        if ((lp = fopen (strconcat(file, ".log"), "w")) == NULL)
        {
                perror("fopen");
                exit(1);
        }
        time(&timer);
        gethostname(hostname, BUFSIZ);

        fprintf(lp, "init   :");
        while (s = *argv++) fprintf(lp, " %s", s);
        fprintf(lp, "\n");
        fprintf(lp, "host   : %s\n", hostname);
        fprintf(lp, "start  : %s", ctime(&timer));
        fprintf(lp, "--------\n");
        fflush(lp);
}

logging::~logging ()
{
        time_t timer;
        time(&timer);
        
        fprintf(lp, "--------\n");
        fprintf(lp, "finish : %s", ctime(&timer));
        fclose(lp);
}

void
init_log ()
{
        if ((lp = fopen (DEFAULT_LOG_FILE, "w")) == NULL)
        {
                perror("fopen");
                exit(1);
        }
}

int
lprintf (const char *fmt, ...)
{
        int n;
        va_list ap;
        
        if (lp == NULL) init_log ();
        va_start(ap, fmt);
        n = vfprintf(lp, fmt, ap);
        va_end(ap);
        fflush(lp);
        return n;
}

中で lprintf() という関数を定義してあるので, これを使うと他の情報もログファイルに自動的に書き込める。 学習の途中で尤度等を記録しておくのによさげです。
もしかするとこういうテクは既に知られているのかも知れないですが, lprintf()みたいなものを絡めたものは多分ないと思うので, 書いてみました。

---

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

-