和が1らしいです。ついにオフィスにもベイズが進出!

という冗談はともかく, 今日はISMの 統計思考院 の 設立シンポジウム でした。 NIIの新井先生の講演「科学の言葉としての数学」が, 教養とビジョン溢れる, 非常に素晴らしいお話でした。
先日の情報・システム研究機構の若手合宿(「クロストーク」＠群馬安中温泉)の 帰りのバスでは, 僕は実は新井先生のすぐ後ろの席に座っていたわけですが, 凄いなあと感動したり。
夜は統計思考院の設立記念式典と祝賀会で, 階段にまでいい匂いがしていました。

石垣島の石垣商工会館で行われた, 情報処理学会第204回自然言語処理研究会から, 先週の23日(祝)に帰ってきました。
NL研は毎回委員会をやると最初に聞いていた気がするのもあって参加したのですが, 実際には委員会はなかったので, 参加必須ではなかったようですが..。 その代わり, 夜に非公式な委員会?がありました。

内容はメモを取りながら聞き, 色々NLPの課題も整理できたので, 良かったような気が します。二日目夜は懇親会の後, 石垣市街のビルの6Fにあるバーで,

　| 高村○ ○船越(本田)
窓| 宮尾○ ○西川(NTT)
　| 持橋○ ○吉川(東工大)

という濃いメンバー(敬称略)で, 泡盛を飲みました。以下写真。(クリックで拡大)

NL研最遠の地。	会場の石垣商工 会館。	ホテルからの風景。 実際の海はもっと トルコブルーです。	西表島, 大原集落の NTT交換局。

先日, 内部で話していてふと思いついた話を実装しようと思い, 式を書いて微分して Newton法を導いてプログラムを書き, 動かしてみたところで, 目的関数が全然凸でない ため初期値依存で全く動かないことが判明。;
MCMCにするしか..と考えたところで, BUGSで計算すればいいのでは, ということに気付き, 久保さん にお聞きして, ついに僕もRと JAGS を使ってみることに。BUGSはWindowsやLinuxでしか動かないので, OSXやSunOS(w で 動かすには, オープンな実装のJAGS等が必要です。

さて, BUGSは多くの統計の場合と同様, {3.28, 4.07, 1.15, -5.92, ..}のような連続値のデータが, ある分布に従ってiidで一度ずつ発生している場合を対象にしているらしく, 自然言語処理で 現れるような, data[i,j] が k回, というような場合は, データを仮想的に

for (k in 1:freq[i,j])
{
    data2[i,j,k] <- data[i,j];
}

などとして増やさなければいけない模様(!)。
このままでは, 頻度が2000のアイテムは計算もメモリも2000回に増えてしまう。 また, BUGSでは一般的な混合モデルが上手く扱えないらしく, 同じ分布の混合分布はそのまま 書けたり, GMMの場合にはパッケージもあるものの, 例えば Gaussian と外れ値 分布としての Cauchy, のような混合分布を扱う方法がない。 *1

途方に暮れたので久保さんにお聞きしてみると, BUGSには尤度を直接書ける Poisson trick と 呼ばれるものがある, とのこと。
これは, パラメータμのポアソン分布から0が出る確率は exp(-μ) であることを利用して,

data[i,j] <- 0; L[i,j] = log(何とか..);
data[i,j] ~ dpois(-L[i,j]);

とすると, 任意の尤度関数が書ける, というもの。 WinBUGS User manual の36ページに載っていて, 一部では有名なトリックらしい。 これを使って, dat[i,j] の頻度が freq[i,j] に書かれている時,

data {
    for (i in 1:N)
    {
        for (j in 1:index[i])
        {
            pseudo[i,j] <- 0;
        }
    }
}
model
{
    ..
    for (i in 1:N)
    {
        ..
        for (j in 1:index[i])
        {
                L[i,j] <- log (q[i] / K + (1 - q[i]) * exp (loggam(dat[i,j] + r[i]) - ..));
                pseudo[i,j] ~ dpois(- L[i,j] * freq[i,j]);
        }
    }
}

などのようにして, 混合分布も, 頻度項も上手く吸収することができた(上の freq[i,j] の使い方を参照)。
NLPでJAGS/BUGSを使う人がもしいれば, 同じ問題に遭遇すると思いますし, 久保さんを少しdisturbしてしまったので, 感謝をこめてここに書いておきます。

---