緑本9章の内容について、RStanでやってみる。
Stanコード
次の内容で作成し、sample.stanという名前で保存した。
data {
int<lower=1> N;
vector[N] x;
int<lower=0> y[N];
real mean_x;
}
parameters {
real beta1;
real beta2;
}
model {
for(i in 1:N){
y[i] ~ poisson(exp(beta1 + beta2 * (x[i] - mean_x)));
}
}コードはPyStanによるはじめてのマルコフ連鎖モンテカルロ法 - Qiitaを参考に少しだけいじった。やはり平均値は引いた方が少し早いようだ。Stanでは変数名にドットを使用できないので注意が必要。
Rコード
library(rstan) 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") # データlistの準備 list.data <- list( x = d$x, y = d$y, N = nrow(d), mean_x = mean(d$x) ) # stanで計算 fit <- stan(file = "sample.stan", data = list.data, iter = 1500, chains = 3, thin = 3, warmup = 100)
緑本メモ - 盆栽日記を参考にした。
結果の確認
> print(fit) Inference for Stan model: sample. 3 chains, each with iter=1500; warmup=100; thin=3; post-warmup draws per chain=467, total post-warmup draws=1401. mean se_mean sd 2.5% 25% 50% 75% 97.5% beta1 1.97 0.00 0.08 1.81 1.92 1.98 2.03 2.13 beta2 0.08 0.00 0.07 -0.06 0.03 0.08 0.13 0.22 lp__ 143.99 0.03 1.02 141.56 143.62 144.31 144.70 144.95 n_eff Rhat beta1 1287 1 beta2 1263 1 lp__ 1315 1 Samples were drawn using NUTS(diag_e) at Mon Aug 15 00:45:10 2016. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).
plot(fit)
結果のプロット
トレースプロットはrstanパッケージのtraceplot()で描ける。バーンイン期間の表示の有無等も設定できるようだ。
traceplot(fit, inc_warmup = TRUE)
また、サンプリング結果はextract()で数値として取り出せるので、これを用いることで事後分布のプロット等も描ける。
post <- extract(fit) old <- par(mfrow = c(1, 2)) plot(density(post$beta1), main = "beta 1") rug(post$beta1) plot(density(post$beta2), main = "beta 2") rug(post$beta2) par(old)
(※2016/08/17追記)
単に事後分布の密度推定をプロットしたいだけであればstan_dens()を使用するのが簡単だった。chain毎の推定結果を重ねて表示するオプションなんかもある。可視化のためのツールは一通り用意されている(?stan_plotを確認するとよい)。
stan_dens(fit1, separate_chains = TRUE)
(※追記ここまで)
テキストの図を再現する
テキストの図9.6については、生態学データ解析 - 本/データ解析のための統計モデリング入門のfig09_06.Rを少し修正すれば描くことができる。
# テキストの図を再現 add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){ lines(d$x, exp(bb1 + bb2 * (d$x - mean(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$beta1)) { add.mean(post$beta1[i], post$beta2[i], lty = 1, col = "#00000004") } points(d$x, d$y) add.mean(median(post$beta1), median(post$beta2), lty = 1, lwd = 2)
