「データ解析のための統計モデリング入門」第9章メモつづき

前回。
rion778.hatenablog.com
これをrstanarmでやってみた。

準備

install.packages("rstanarm")

割と時間かかった。

library(rstan)
library(rstanarm)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# データダウンロード
download.file(destfile = "ch9.Rdata", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData")
load("ch9.Rdata")

実行

fit <- stan_glm(y ~ x, family = poisson(), data = d)

細かい設定は省略したけど1行で済むのは楽。あと早い。

プロット

rstanに入っている作図関数がそのまま使える。

plot(fit)
stan_dens(fit, separate_chains = TRUE) + geom_rug()
stan_trace(fit)

出力は特に違いがないので略。
stan_traceにinclude_warmup=TRUEを指定した時の挙動が若干おかしい気はする(thinが無視されている様子?)。

出力

デフォルトの出力はかなり要約されている。summaryで少し詳しい物が見られるが、$stanfitの中にrstanの時と同じ形式のやつが入っていたので、例えばextract()を使いたい時なんかはこれを使えば良い様だ。

fit                # デフォルト
summary(fit)       # 少し詳しく
fit$stanfit        # rstanと同じで

テキストの図を再現する

上述のように、$stanfitをextract()してやれば良い。変数名が何故かalphaとbetaだ。

post <- extract(fit$stanfit)
add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){
  lines(d$x, exp(bb1 + bb2 * d$x), lty = lty, lwd = lwd, ...)
}
plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
for (i in 1:length(post$alpha)) {
  add.mean(post$alpha[i], post$beta[i], lty = 1, col = "#00000004")
}
points(d$x, d$y)
add.mean(median(post$alpha), median(post$beta), lty = 1, lwd = 2)

これも出力に差は無いので略。