読者です 読者をやめる 読者になる 読者になる

Knowledge As Practice

JAIST(東京)でサービス経営の研究をしている社会人大学院生の研究・勉強メモ(統計分析多め)。

統計学エンドユーザーのための JAGS 入門

ベイズデータ分析 統計分析 ノンプログラマーのためのR入門

岩波データサイエンス vol.1 のどこかのページに「初心者はJAGS(BUGS)」と書いてあったので*1、こつこつ学び始めています。世の中は Stan 大人気ですが、気にしない…。まず JAGS に慣れてから Stan をやります…。本エントリーは、前回のエントリーと多くがかぶっています。

 
今回は、次のページを参考にさせていただいています。

生態学のデータ解析 - R で JAGS (rjags 編)
http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html

生態学のデータ解析 - JAGS
http://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html

 
分析に使用したのは、JAGS 4.0.0 および rjags 4.4 です。

 

JAGS と rjags のインストール

JAGS と rjags はそれぞれ別にインストールします。JAGS 3.0 系 をインストールしていて、rjags 4.4 にしていた場合は動きません。JAGS も4にする必要があります。

 
JAGS はこちらから .exe をダウンロードできます。
http://sourceforge.net/projects/mcmc-jags/files/

rjags はいつもの install.packages でOKです。

 

分析の手順

1.データの準備

今回は先ほどの「RでJAGS(rjags 編)」のデータを使用しました。上記ページにある csv ファイルから、ページのとおりにデータをセットします。JAGSに渡すデータはリスト形式である必要があります。

 

2.モデルを作る

「RでJAGS(rjags 編)」にあるモデルを自分なりにいじってみました。いまのJAGSマニュアルによると、行末にセミコロンは不要のようです。モデルを作る前に図解しました。これも岩波データサイエンスにそうしろ、と書いてあったように思います。分布の形はテキトーです*2

f:id:hikaru1122:20151120002559j:plain

 
さて、モデルです。ロジスティック回帰でリンク関数はロジットリンクです。詳しくはみどりぼん第6章を参照してください。

Model <-"
model{
    for (i in 1:N) {
        Positive[i] ~ dbin(p[i], Total[i])
        logit(p[i]) <- z[i]
        z[i] <- beta[1] + beta[2] * Logdose[i]
    }
    beta[1] ~ dnorm(0, tau)
    beta[2] ~ dnorm(0, tau)
    tau ~ dgamma(1.0E-6, 1.0E-6)
}
"

 

追記2015年11月20日
そういえば、JAGS の dnorm の分散は逆数でした。
tau <- 1/(sd*sd)
sd ~ dgamma(1.0E-3, 1.0E-3)
がより正しい書き方なのかな…?
また、twitter にてアドバイスをいただきました。
一様分布や半コーシー分布のほうがよい、とのことです。
ありがとうございます。

 
いちばんの変更点は tau の事前分布です。一度、同じように dgamma(1.0E-3, 1.0E-3) でやってみたのですが、うまく行きませんでした。代わりに、dgamma(1.0E-6, 1.0E-6) を使いました。

 
その理由は2つあります。
(1)JAGS 4.0 マニュアルに載っていた例で dgamma(1.0E-6, 1.0E-6) が使われていたこと
(2)前回のエントリーで dgamma(1.0E-6, 1.0E-6) を使ってやり直したら、Stan とほぼ同じ結果が得られたこと
です。1.0E-3 と 1.0E-6 で結果が変わるのってシビア。事前分布って難しい。

 
モデルを書いたら、書き出します。

writeLines(Model, con ="testModel.txt") # 拡張子は特に不要ですが、いちおう。

 

3.MCMCする

モデルを jags オブジェクトにしていきます。

# rjags のロード。
library(rjags)
# jags モデルを作る。
testModel <- jags.model(file = "testModel.txt", data = list.data, n.chain = 3)

 
jags.model の file = は別ファイルである必要があります。なので、先ほど書きだした次第です。次にバーンイン、そして本格的なサンプリングを行います。n.chain は2以上でないと、後で Gelman & Rubin 収束診断指標のグラフが描けないようなので、3にしています。初期値は設定していません。まだ理解できていないので。

 
では、本格的なMCMCサンプリングを行います。

# バーンイン期間のための試運転。
update(testModel, n.iter = 1000)
# 本番。MCMC サンプルを kekka に収める。
kekka <- coda.samples(
            testModel,
            variable.names = c("beta[1]", "beta[2]", "tau"),
            n.iter = 20000)

 
coda.samples の variable.names では結果で見たいパラメータを文字列として指定します。これで、下ごしらえは終わりました。rjags パッケージで使うのは、jags.model, update, coda.samples の3つだけです。writeLines は普通に使えます。

 

結果を見る。

収束診断をします。codamenu() でメニューを選びながら、確認します。引数はいらないです。メニューで「2→kekka(MCMCサンプリングを格納したオブジェクト)→2→2」とするとGelman & Rubin の指標が見られます。1.1 より下なので、収束しているようです。他にもcodamenu()でいろいろ指標を選べます。

Potential scale reduction factors:
        Point est. Upper C.I.
beta[1]       1.01       1.02
beta[2]       1.01       1.02
tau           1.00       1.00

 
図でも確認します。plot(kekka)です。

f:id:hikaru1122:20151120011135p:plain

 
推定値は summary(kekka) で確認します。

> summary(kekka) 
            Mean      SD  Naive SE Time-series SE
beta[1] -5.99179 1.10916 0.0045281       0.053419
beta[2]  1.18754 0.19647 0.0008021       0.009398
tau      0.06042 0.07121 0.0002907       0.001165

 
「RでJAGS(rjags 編)」ページに書いてある「真の値」は、 beta[1] = -5.0、beta[2] = 1.0 です。けっこう近づきました。tau だけはダメでした。図示します。

f:id:hikaru1122:20151120011149p:plain

 
少しずつ、JAGS に慣れてきたようです。今後、みどりぼんや岩波データサイエンスの例にも取り組んでいきたいと思います。

 

参考文献

JAGS マニュアル http://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/

rjags マニュアル https://cran.r-project.org/web/packages/rjags/rjags.pdf

生態学のデータ解析 - R で JAGS (rjags 編) http://hosho.ees.hokudai.ac.jp/~kubo/ce/RtoJags.html

生態学のデータ解析 - JAGShttp://hosho.ees.hokudai.ac.jp/~kubo/ce/JagsMisc.html

岩波データサイエンス

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

みどりぼん

DBDA2

Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan

Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan

*1:出張中でいま手元になく、どのページかわかりません…。

*2:矢印は上向きのほうがいいのかな。

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 - 継承 4.0 国際 ライセンスの下に提供されています。