前回。
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)
これも出力に差は無いので略。