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

Knowledge As Practice

JAIST(東京)で Transformative Service Research に取り組んでる社会人大学院生の研究・勉強メモ

JAGS と Stan と lm

ここ最近、GLM や GLMM の学び直しをしていました。みどりぼんも2周目になり、6~7割は理解できるようになったと思います。そこで、唐突ですが3ヶ月ほど前に Stan&Rstan でやった、単回帰分析を JAGS&rjags でやります。元ネタ(?)はこちら。 hikaru1122.hatenadiary.jp

 
JAGS を使うにあたって参考にした書籍はこちら(以下、DBDA2 と呼びます)。また、JAGSのユーザーマニュアルもざっと読みました。

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

sourceforge.net

 
DBDA2 の第8章でやっと JAGS を使った分析例が出てきます。DBDA2 による JAGS を使ったベイズ推定の手順は次のとおりです*1。勝手に3つの段階「準備する→MCMCする*2→結果を見る」にわけてみました。

 

<準備する>

1.分析するデータを準備する。
2.リスト形式でまとめる。…A
3.モデルを作る。
4.モデルを書き出す。…B
5.(初期値を設定する…C)
※5は必ずしも必要ではない。

 

MCMCする>

1.JAGSにすべての情報(A,B,C)を渡す。
→ rjags::jags.model を使う。

2.バーンイン期間のために空回しする。 → rjags::update を使う。

3.結果を記録する → rjags::coda.samples

 

<結果を見る>

・収束しているかどうか確認する。
→ 著者が用意してくれた diagMCMC を使う*3

・事後分布を確認する。
→ mode, median, mean を確認。mode がいちばんいいけど、不安定なことももある。median がたいていいちばん安定している。
→ 著者が用意してくれた plotMCMC を使う*4

 
では、やってみます。分析するデータは「アイスクリーム統計学」からお借りします。 http://kogolab.chillout.jp/elearn/icecream/chap4/sec1.html

 

<準備する>
icecream <- read.delim("icecream.txt")
kyaku <- icecream$kyaku
kion <- icecream$kion
Ntotal <- length(kyaku)

# JAGS に渡すデータをリスト形式で作る。
bunsekiData <- list(
    kyaku = kyaku,
    kion = kion,
    Ntotal = Ntotal
)

# モデルを作る。
Model <-"
model{
    for(i in 1:Ntotal){
    kyaku[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta * kion[i]
    }
    alpha ~ dnorm(0, 1.0E-4) #事前分布はこれでいいの?
    beta ~ dnorm(0, 1.0E-4) #事前分布はこれでいいの?
    tau ~ dgamma(1.0E-3, 1.0E-3) #事前分布はこれでいいの?
}
"

#モデルを書き出す。
writeLines(Model, con = "icecreamModel.txt")

 

MCMCする>
#rjags をロード。
library(rjags)

# jagsオブジェクトを作る。初期値は設定せず、おまかせ。
haahaa <- jags.model(file = "icecreamModel.txt",
    data = bunsekiData, n.chains = 3, n.adapt = 500)

# バーンイン期間のための空回し。
update(haahaa, n.iter = 5000)

# MCMCサンプルを kekka に収める。知りたいパラメータは alpha(切片)と beta(偏回帰係数)。
kekka <- coda.samples(haahaa, variable.names = c("alpha", "beta"), n.iter = 100000)

 

<結果を見る>
source("DBDA2E-utilities.R") # これは DBDA2 の付録で付いてくる。

# DBDA2 著者作成の関数で収束診断。
diagMCMC(codaObject = kekka, parName = "alpha")
diagMCMC(codaObject = kekka, parName = "beta")

# DBDA2 著者作成の関数で事後分布を見る。mean, mode, median をそれぞれ見てみる。
par(mfrow=c(2,3))
plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "mean")
plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "mode")
plotPost(kekka[,"alpha"], main = "alpha", xlab = bquote(alpha), cenTend = "median")
plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "mean")
plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "mode")
plotPost(kekka[,"beta"], main = "beta", xlab = bquote(beta), cenTend = "median")

 
収束はしているようです(下図のそれぞれ左上と左下)。 f:id:hikaru1122:20151118171650p:plain f:id:hikaru1122:20151118171705p:plain

 
切片と係数は次のとおりです。それぞれ3つ出してみます。 f:id:hikaru1122:20151118171852p:plain

 
median を採用するとしたら、回帰式は

{ \displaystyle
客数 = -145 + 14.6*気温
}

ということですよね。

 
同じデータを lm でやった結果は次のとおりです。

> kekka_lm <- lm(kyaku ~ kion, data=icecream)
> summary(kekka_lm)

Call:
lm(formula = kyaku ~ kion, data = icecream)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.969 -17.709  -1.218  17.413  51.031 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -229.98      73.79  -3.117  0.00596 ** 
kion           17.25       2.30   7.499 6.08e-07 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 29.54 on 18 degrees of freedom
Multiple R-squared:  0.7575,    Adjusted R-squared:  0.744 
F-statistic: 56.23 on 1 and 18 DF,  p-value: 6.082e-07

Stan での結果も再掲します。

         mean se_mean    sd    2.5%     25%     50%     75%  97.5% n_eff Rhat
alpha -230.75    0.57 80.23 -389.66 -282.82 -230.41 -178.89 -72.06 19789    1
beta   17.27    0.02  2.50   12.32   15.66   17.26   18.90  22.23 19771    1

 
Stan&Rstan でやると lm とほぼ同じ結果が得られたんですが、今回の JAGS だとなんか違います。事前分布の設定が悪かったのかな…。それとも初期値の設定が必要だったのか。それとも分析するデータが少なすぎるとか。やっとモデルを作る、見るのに慣れてきたと思ったら、また次の難問が出てきました。厳しい…。

 
しかし、Stan よりは JAGS のほうがとっつきやすい気がします。しばらく慣れ親しんでみようと思います。

*1:Stanでも大きくは変わらないと思うのですが、念のため

*2:MCMCするという謎な日本語は許してください。

*3:ダウンロード先は、https://sites.google.com/site/doingbayesiandataanalysis/software-installation

*4:上に同じ。

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