Created
May 24, 2019 18:17
-
-
Save Gedevan-Aleksizde/88cc7c3518128471d67e99dbc98bac1c to your computer and use it in GitHub Desktop.
Tokyo.R 78'
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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