JAGSを使ってギブスサンプリングを試してみた (1)

  • ????????????????????

前から気になってた JAGS を試してみました。Just Another Gibbs Sampler を略して JAGS。Bugs ファミリーの中では現在、最もアクティブに開発がすすめられているようです。

ここでは、とにかく何も知らない状態から R+JAGS で簡単なモデル(ガウスノイズな線形回帰)の推定までできるように最低限のことをメモしました。

準備

JAGS本体のインストール

JAGS は R とは分離された独立なソフトウェアなので R の CRAN 経由で一発インストールというわけにはいきません。ですが、Ubuntu の場合は apt で一発です。他のOSは試していないので、インストールしていない人は別途調べてください。

sudo apt-get install jags

rjagsのインストール

JAGSをRから呼び出すためには rjags というパッケージをインストールする必要があります。apt 使うなら

sudo apt-get install r-cran-rjags

もしくはRのプロンプトで

install.pacakges("rjags")

とすればOK。

テスト

テスト問題として

\begin{align*}
y_i&\sim \text{Normal}(\mu_i,\sigma^2)\\
\mu_i&=\alpha+\beta x_i
\end{align*}
というデータの発生メカニズムを仮定します。これに対応した次のような「モデルファイル」をつくって、適当に名前をつけておきます(lm.jagsとする)。

model {
  for(i in 1:N){
    mu[i] <- alpha + beta*x[i]
    y[i] ~ dnorm(mu[i],1/(tau^2))
  }
  alpha ~ dnorm(0,0.0001)
  beta ~ dnorm(0,0.0001)
  tau ~ dgamma(0.1,0.1)
}

見た感じ、Rのコードにかなり似ているけど、実際には細かい点で異なります。dnormが正規分布なのは同じだけど、第一引数は存在しないし、JAGSでは分散ではなくて精度(precision;分散の逆数)を指定することになっているので注意。Rを知っていればなんとなくわかってしまう構文は素晴らしいですが、上で挙げたような注意の他に、分布の名前自体がちがっていたりすることもあるのでハマる可能性もあります。ちなみに、for ブロックの後ろの3行は未知パラメータの事前分布を設定しています。

データの生成

\(\alpha=2, \beta=3\) としてデータを発生させます。

# determine the number of samples
N <- 50
# data generation
x <- rnorm(N) 
y <- 2 + 3*x + rnorm(N,0.1)

JAGSをRから呼び出す

モデルファイルを保存したら、Rで以下のスクリプトを走らせます。すると下のようなグラフが表示されるはずです。coda[1]でプロットさせているのでMCMCglmmなど他のRのベイズパッケージを使った時と同じ図になります。

library(rjags)
library(coda)

# compiling and initializing
lm.mcmc <- jags.model("lm.jags",
                       data=list('x'=x, 'y'=y, 'N'=N),
                       n.chains=4,  # the number of parallel chains
                       n.adapt=100)

# mcmc sampling
update(lm.mcmc,1000) # burn-in
posterior <- coda.samples(lm.mcmc, c('alpha','beta'), 2000)
# plotting the result
plot(posterior) # plot.coda is called


こんな感じで推定出来ました。すこしスクリプトについて書いておくと

  • jags.model 関数の第一引数では先程のモデルファイルを指定する。
  • 第二引数 data では JAGS に渡したいRオブジェクトをリスト形式で与える。クオーテーションする必要があるので注意。
  • n.chains で並列パスの数を指定
  • update で burn-in している
  • coda.samples で指定した数のサンプルを coda 形式で取得する。

まとめといろいろと蛇足

階層ベイズでも何でもないつまらない例を使って、ほんとに最低限のJAGSの使い方について書きました。たぶん続編があるはずです。次は混合ポアソン回帰でも試してみようかな。

わりとどうでもいいことですが、’J’AGS という名前(←JRubyみたいな意味で。)とクロスプラットフォームを謳っていることから Java で書かれていると思っていましたが、c++で実装されているようです。この勘違いのせいでずいぶん長い間敬遠してきた、という経緯があります。

ちなみに、これまで自分でMCMC実装するときは

  • cで書く
  • Rでものすごく遅いやつを書く(ちかいうちにこのブログにも登場予定)
  • MCMCglmmみたいなRのパッケージを使う

という選択肢の中から選んでいたのですが、これにJAGSを追加して比較するとこんな感じでしょうか?(いろいろと異論はあるかもしれませんが)

自由度 実行スピード 実装の手軽さ
c言語 ×
R ×
R-packages × ○ (?)
JAGS ○ (?)

おしまい。

はてなブックマーク - JAGSを使ってギブスサンプリングを試してみた (1)
Pocket

  1. [1] RでMCMCの結果を可視化するときのほぼ事実上の標準形式