Tokyo.R 78'

bayesplot_tutorial.R

require(rstan)
require(bayesplot)
require(ggplot2)
require(ggthemes)
require(patchwork)
require(tidyverse)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

#### rstan でモデルをコンパイルしサンプリング ####

source_8school <- '
// https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started-(Japanese)
data {
  int<lower=0> J;         // 学校の数
  real y[J];              // 推定されている教育の効果
  real<lower=0> sigma[J]; // 教育の効果の標準誤差
}

parameters {
  real mu;                // 処置の効果（全体平均）
  real<lower=0> tau;      // 処置の効果の標準偏差
  vector[J] eta;          // 学校ごとのスケール前のバラつき
}

transformed parameters {
  vector[J] theta = mu + tau * eta;        // 学校ごとの処置の効果
}

model {
  target += normal_lpdf(eta | 0, 1);       // 事前分布の対数密度
  target += normal_lpdf(y | theta, sigma); // 対数尤度
}
'

model_8school <- stan_model(model_code=source_8school)

# 公式の例よりもデータとパラメータを増やす
set.seed(42)
schools_dat <- list(J=100)
schools_dat[['sigma']] <- with(schools_dat, rchisq(n=J, df=2))
schools_dat[['y']] <- with(schools_dat, rnorm(n=J, sd=sigma))

samples <- sampling(
  model_8school,
  data=schools_dat,
  chains=4,
  iter=2000, warmup=500
  )

#### 事後診断 ####

print(samples)
samples_as_array <- as.array(samples)
summary_samples <- summary(samples)$summary

# Rhat の確認
mcmc_rhat(bayesplot::rhat(samples))


# 有効サンプルサイズの確認
mcmc_neff(bayesplot::neff_ratio(samples))

# トレースプロット
mcmc_trace(samples_as_array, pars='mu')
mcmc_trace(samples_as_array, pars='tau')

# 連番のパラメータを分割表示するのは少し面倒
# 正規表現も使えるが, 今回は不適当
mcmc_trace(samples_as_array, regex_pars='^eta')

# こうするのが楽?
for(x in with(schools_dat, split(1:J, ceiling(seq_along(1:J)/20)))){
  print(mcmc_trace(samples_as_array, pars=paste0('eta[', x, ']')))
}
for(x in with(schools_dat, split(1:J, ceiling(seq_along(1:J)/20)))){
  print(mcmc_trace(samples_as_array, pars=paste0('theta[', x, ']')))
}

# 自己相関係数
mcmc_acf_bar(samples_as_array, pars='mu')
mcmc_acf_bar(samples_as_array, pars='tau')

for(x in with(schools_dat, split(1:J, ceiling(seq_along(1:J)/20)))){
  print(mcmc_acf_bar(samples_as_array, pars=paste0('eta[', x, ']')))
}
for(x in with(schools_dat, split(1:J, ceiling(seq_along(1:J)/20)))){
  print(mcmc_acf_bar(samples_as_array, pars=paste0('theta[', x, ']')))
}


mcmc_nuts_divergence(bayesplot::nuts_params(samples), bayesplot::log_posterior(samples))

#### 収束しないダメなモデル ####

# パラメータが一意にならない
# https://rpubs.com/piroyoung/stan01
model1 <- '
parameters{
  real a;
  real b;
}
model{
  a ~ normal(0,10000);
  b ~ normal(0,10000);

  0.0 ~ normal(a + b,1.0);
}
'
fit <- stan(model_code=model1)
mcmc_trace(as.array(fit), pars=c('a', 'b'), facet_args=list(ncol=1), size=1) +
  theme(text=element_text(size=20))

mcmc_pairs(as.array(fit), pars=c('a', 'b')) +
  theme(text=element_text(size=20))

# label switching
model1 <- '
data {
  int N ;
  vector[N] y ;
}
parameters{
  real<lower=0, upper=1> theta ;
  real mu[2] ;
  real<lower=0> sigma[2] ;
}
model{
  mu ~ normal(0, .6) ;
  sigma ~ uniform(0, 1) ;
  target += log_mix(theta, normal_lpdf(y | mu[1], sigma[1]), normal_lpdf(y | mu[2], sigma[2])) ;
}
'

set.seed(42)
data_labswich <- list(N=100)
data_labswich[['y']] <- with(data_labswich, 0.3 * rnorm(N, -4, 1) + 0.4 * rnorm(N, 0, 2))

fit <- stan(model_code=model1, data=data_labswich, chains=2)

print(fit)
mcmc_trace(as.array(fit), regex_pars='mu', facet_args=list(ncol=1), size=1) +
  theme(text=element_text(size=20))

#### スライド用の画像作成 ####

mcmc_rhat(bayesplot::rhat(samples)[1:20], size = 5) + 
  theme(text = element_text(size=20))
mcmc_neff(bayesplot::neff_ratio(samples)[1:20], size=5) + 
  theme(text = element_text(size=20))
mcmc_trace(samples_as_array, pars=paste0('eta[', 1:8, ']')) +
  theme(text=element_text(size=20))
mcmc_acf_bar(samples_as_array, pars='tau') +
  theme(text=element_text(size=20))

# あれ
require(png)
# ggplot オブジェクトと同じように扱えるので, scale_color_*() で色変更もできる.
# patchwork で画像をつなげることも可能
g1 <- mcmc_trace(samples_as_array, pars='mu') + labs(y=expression(mu), title='Trace plot') +
  theme(axis.title.y=element_text(angle=0, vjust=.5)) +
  scale_color_brewer(palette='GnBu', direction=-1)
catap <- ggplot(data.frame(x=1:100, y=1:100), aes(x=x, y=y)) +
  annotation_raster(raster=readPNG('hairycatepiller.png'), -Inf, Inf, -Inf, Inf) +
  labs(title='Hairy caterpillar (from Costa Rica)') + theme(axis.title=element_blank())
g1 / catap
ggsave('traceplot_and_caterpillar.pdf', g1 / catap)