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

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

緑本9章の内容について、RStanでやってみる。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

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)

f:id:Rion778:20160815014841p:plain

結果のプロット

トレースプロットはrstanパッケージのtraceplot()で描ける。バーンイン期間の表示の有無等も設定できるようだ。

traceplot(fit, inc_warmup = TRUE)

f:id:Rion778:20160815015315p:plain
また、サンプリング結果は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)

f:id:Rion778:20160815020229p:plain
(※2016/08/17追記)
単に事後分布の密度推定をプロットしたいだけであればstan_dens()を使用するのが簡単だった。chain毎の推定結果を重ねて表示するオプションなんかもある。可視化のためのツールは一通り用意されている(?stan_plotを確認するとよい)。

stan_dens(fit1, separate_chains  = TRUE)

f:id:Rion778:20160817204956p:plain
(※追記ここまで)

テキストの図を再現する

テキストの図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)

f:id:Rion778:20160815020806p:plain