Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Tokyo.R 78'
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.