Introduction to bayesplot (mcmc_ series)

bayesplot packageについて
使う前の準備
Stan model
Stanを回す
bayesplotを使う [本番]
関数を探したいときに便利な関数
とりあえずこのパッケージで遊びたい
まとめ
今回使用した全Rコード

こちらは、Stan advent (Link) 2日目の記事になります。
お前誰や？？→→HPに行く

昨日の小杉先生M1の記事は面白かったですね!!!
今日は、baysplotの紹介をしたいと思います。

bayesplot packageについて

Stan楽しいですね(*´Д`)ﾊｧﾊｧ

ベイズ推定をするのも楽しいですが、ただベイズ推定して終わりではなく、結果の解釈が一番大事であり、その前にやるべきこともあります。

Stanでパラメータ推定を行う場合MCMC(HMC)を行うと思いますが、推定がきちんと終わっているか(収束しているか)を確認することは非常に大事です。

推定の確認方法として以下をみることが一般的かと思います。
- パラメータ推定が収束しているかを確認する(Gelman-Rubin診断の\(\hat{R}\)等)
- パラメータ推定(MCMC)が安定しているか確認する(effective sample size(n_eff/ESS)やauto correlation等※)
- 事後分布を確認
- 事後予測チェックをする

※ ESSの計算にauto correlationを使っているので、別物ではありませんが…
(そういや最近、Stanのauto correlationの近似方法が変わりましたね。。。)

Stanには、既にggplot2ライクな事後処理用の関数が標準装備されています。
ここら辺に関しては、以下の資料をご覧ください。

Stan便利な事後処理関数SlideShare

これとは別にStan公式が出しているbayesplotパッケージというものがあって、これも同じようなことができます。
既に日本語文献として小杉先生のサイトで紹介されていました(こちら)。
上記も併せてご覧ください。この記事を書くうえでも参考にしています。

bayesplotパッケージでは、大きく二つの関数系に分かれています
1. mcmc_ シリーズ(収束や事後分布を確認するための関数群)
2. ppc_ シリーズ(事後予測チェックを行うのに便利な関数群)

この記事では、1.のmcmc_シリーズの説明を行います。
2.については誰かもしくは未来の自分がやってくれることでしょうｗ

使う前の準備

とりあえず適当にStanを回します。
今回は、デモで最もよく使われるであろう8schoolモデルを回します。
今回は特に事後処理の説明のみするので、ここではモデルがどうなっているかを見る必要はありません。

# setwd('')

rm(list = ls())
library(tidyverse)  ## ほぼパイプのために読み込んでいます。
library(rstan)
library(bayesplot)  ## こいつが今回の主役です。
rstan_options(auto_write = TRUE)  ## なくても構いません
options(mc.cores = parallel::detectCores())  ## なくても構いません
options(max.print = 99999)  ## なくても構いません

そして、使うデータを準備します。

schools_dat <- list(J = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 
    10, 16, 11, 9, 11, 10, 18))

Stan model

今回回すStanモデルはこちらです。

data {
  int<lower=0> J;
  real y[J];
  real<lower=0> sigma[J];
}

parameters {
  real mu;
  real<lower=0> tau;
  real eta[J];
}

transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
}

model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}

上記を自分で動かすなら、model.stanと保存して、R上で以下を実行してください。
モデルをダウンロードする

model <- stan_model("model.stan")

Stanを回す

以下でStanが回ります。

fit <- sampling(object = model, data = schools_dat, iter = 1000, chains = 4)

## 結果 ##
fit

Inference for Stan model: 4a38cae9bad394789360bce03f17c3ed.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

           mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu         8.05    0.14 4.96  -1.50   4.82   8.14  11.14  18.17  1243 1.00
tau        6.25    0.16 5.19   0.25   2.32   5.17   8.85  19.03  1042 1.00
eta[1]     0.40    0.02 0.94  -1.54  -0.21   0.40   1.05   2.18  1824 1.00
eta[2]     0.00    0.02 0.94  -1.86  -0.62  -0.01   0.62   1.90  1555 1.00
eta[3]    -0.20    0.02 0.97  -2.14  -0.84  -0.19   0.43   1.74  1794 1.00
eta[4]    -0.06    0.02 0.87  -1.82  -0.61  -0.05   0.49   1.65  1968 1.00
eta[5]    -0.36    0.02 0.90  -2.14  -0.95  -0.40   0.21   1.50  1719 1.00
eta[6]    -0.23    0.02 0.91  -1.98  -0.83  -0.26   0.36   1.65  1983 1.00
eta[7]     0.32    0.02 0.86  -1.45  -0.24   0.36   0.89   1.97  1958 1.00
eta[8]     0.07    0.02 0.94  -1.78  -0.58   0.09   0.70   1.94  1839 1.00
theta[1]  11.54    0.21 8.35  -1.63   6.19  10.36  15.45  31.23  1550 1.00
theta[2]   7.93    0.14 6.47  -5.21   3.95   7.76  11.95  20.60  2051 1.00
theta[3]   6.25    0.18 7.90 -12.11   2.25   6.67  11.21  20.69  1830 1.00
theta[4]   7.60    0.14 6.36  -6.39   3.84   7.86  11.52  19.67  2058 1.00
theta[5]   5.38    0.14 6.43  -8.40   1.73   5.78   9.52  16.85  2041 1.00
theta[6]   6.23    0.15 6.63  -8.51   2.50   6.66  10.51  18.56  1995 1.00
theta[7]  10.67    0.15 6.70  -0.99   6.28  10.21  14.35  26.00  2031 1.00
theta[8]   8.85    0.18 7.69  -5.32   4.19   8.50  13.01  24.97  1833 1.00
lp__     -39.73    0.10 2.65 -45.60 -41.39 -39.43 -37.81 -35.34   728 1.01

Samples were drawn using NUTS(diag_e) at Sat Dec 01 12:41:12 2018.
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).

bayesplotを使う [本番]

\(\hat{R}\)の確認

一般的に\(\hat{R}<1.10\)で収束とみなすと思います。

## デフォルト not bayesplot
stan_rhat(fit)

stan_rhat(fit)$data  ## dataとしてとってきたい

              stat
mu       0.9992001
tau      1.0038763
eta[1]   1.0004295
eta[2]   0.9995890
eta[3]   1.0003983
eta[4]   0.9994688
eta[5]   1.0001802
eta[6]   0.9996924
eta[7]   1.0003880
eta[8]   0.9997078
theta[1] 0.9995893
theta[2] 0.9989432
theta[3] 1.0008485
theta[4] 0.9996523
theta[5] 1.0009598
theta[6] 0.9997703
theta[7] 0.9994518
theta[8] 1.0007286
lp__     1.0056553

## bayesplot

## rhatがdataとして出力される
rhat(fit)

       mu       tau    eta[1]    eta[2]    eta[3]    eta[4]    eta[5] 
0.9992001 1.0038763 1.0004295 0.9995890 1.0003983 0.9994688 1.0001802 
   eta[6]    eta[7]    eta[8]  theta[1]  theta[2]  theta[3]  theta[4] 
0.9996924 1.0003880 0.9997078 0.9995893 0.9989432 1.0008485 0.9996523 
 theta[5]  theta[6]  theta[7]  theta[8]      lp__ 
1.0009598 0.9997703 0.9994518 1.0007286 1.0056553

## rhatが図として出力
mcmc_rhat(rhat(fit))

## 表示した図に関するdata
mcmc_rhat_data(rhat(fit))

# A tibble: 19 x 5
   diagnostic parameter value rating description 
   <chr>      <fct>     <dbl> <fct>  <chr>       
 1 rhat       theta[2]  0.999 low    hat(R) > 1.1
 2 rhat       mu        0.999 low    hat(R) > 1.1
 3 rhat       theta[7]  0.999 low    hat(R) > 1.1
 4 rhat       eta[4]    0.999 low    hat(R) > 1.1
 5 rhat       eta[2]    1.000 low    hat(R) > 1.1
 6 rhat       theta[1]  1.000 low    hat(R) > 1.1
 7 rhat       theta[4]  1.000 low    hat(R) > 1.1
 8 rhat       eta[6]    1.000 low    hat(R) > 1.1
 9 rhat       eta[8]    1.000 low    hat(R) > 1.1
10 rhat       theta[6]  1.000 low    hat(R) > 1.1
11 rhat       eta[5]    1.00  low    hat(R) > 1.1
12 rhat       eta[7]    1.00  low    hat(R) > 1.1
13 rhat       eta[3]    1.00  low    hat(R) > 1.1
14 rhat       eta[1]    1.00  low    hat(R) > 1.1
15 rhat       theta[8]  1.00  low    hat(R) > 1.1
16 rhat       theta[3]  1.00  low    hat(R) > 1.1
17 rhat       theta[5]  1.00  low    hat(R) > 1.1
18 rhat       tau       1.00  low    hat(R) > 1.1
19 rhat       lp__      1.01  low    hat(R) > 1.1

## histgramとして表示(stan_rhatに近い形で出力)
mcmc_rhat_hist(rhat(fit))

## 特定のパラメータのR hat をみたい。
stan_rhat(fit, pars = c("mu", "eta"))  ## stanの関数

mcmc_rhat(rhat(fit, pars = c("mu", "eta")))  ## bayesplotの関数

使う上でのポイントは,パラメータを指定するしないにかかわらず,rhat()関数にいれてから、mcmc_rhat()にいれて出力
~~stan_rhat()でええやん~~

traceplotを見る

bayesplotを使うためには、stanの結果ファイルを配列(array)型にする必要があります。
stanオブジェクトを配列型にすると、[iter,chains,pars]の3次元配列になります。

## デフォルト
stan_trace(fit, inc_warmup = T)  ## inc_warmupでwarmup区間指定可能

traceplot(fit, inc_warmup = T)

## bayesplot
samples <- as.array(fit)
mcmc_trace(samples, pars = "mu")  ## n_warmupでwarmup区間指定可能

## 特定のchainを特に見たい(目立たせて確認したい)
mcmc_trace_highlight(samples, pars = "mu", highlight = 2)

## 事後分布と同時に確認したい(パラメータ多いと重いので指定推奨)
mcmc_combo(samples, pars = c("mu"))

ポイントは、stanオブジェクトを配列型に変換してからわたすこと。
また、n_warmupの挙動が若干謎。おそらく、warmupとそのあとを同時にプロットできなそう？

~~stan_trace()でええやん~~

auto correlation(自己相関)を確認

自己相関を見るためにも、配列型にする必要があります。先ほどのas.array(fit)をやるのを忘れずに。

## デフォルト
stan_ac(fit, pars = c("mu", "tau"))

## chain別
stan_ac(fit, pars = c("mu", "tau"), separate_chains = T)

## bayesplot

## 折れ線表記
mcmc_acf(samples, pars = c("mu", "tau"))

## 棒グラフ表記
mcmc_acf_bar(samples, pars = c("mu", "tau"))

~~stan_ac()でええやん~~

ESS(n_eff)を確認

続いて、実行(有効)サンプルサイズを確認したいと思います。
今回は、stanオブジェクトをそのまま使うことができます。
一般的にこの値が10％(0.1)以下になると自己相関が高くMCMCがうまくいっていないと判断されます。

## デフォルト
stan_ess(fit)

stan_ess(fit)$data  ## dataとして出力

              stat
mu       0.6212756
tau      0.5211331
eta[1]   0.9120741
eta[2]   0.7774456
eta[3]   0.8970617
eta[4]   0.9837958
eta[5]   0.8594550
eta[6]   0.9912680
eta[7]   0.9790504
eta[8]   0.9192841
theta[1] 0.7748646
theta[2] 1.0252659
theta[3] 0.9149762
theta[4] 1.0290931
theta[5] 1.0205462
theta[6] 0.9975286
theta[7] 1.0152921
theta[8] 0.9165216
lp__     0.3641157

## bayesplot

## n_effの値をdataとして確認
neff_ratio(fit)

       mu       tau    eta[1]    eta[2]    eta[3]    eta[4]    eta[5] 
0.6212756 0.5211331 0.9120741 0.7774456 0.8970617 0.9837958 0.8594550 
   eta[6]    eta[7]    eta[8]  theta[1]  theta[2]  theta[3]  theta[4] 
0.9912680 0.9790504 0.9192841 0.7748646 1.0252659 0.9149762 1.0290931 
 theta[5]  theta[6]  theta[7]  theta[8]      lp__ 
1.0205462 0.9975286 1.0152921 0.9165216 0.3641157

## 図で確認
mcmc_neff(neff_ratio(fit))

mcmc_neff_data(neff_ratio(fit))  ## 図で出力したdataを出力

# A tibble: 19 x 5
   diagnostic parameter value rating description    
   <chr>      <fct>     <dbl> <fct>  <chr>          
 1 neff_ratio lp__      0.364 ok     N[eff]/N <= 0.5
 2 neff_ratio tau       0.521 high   N[eff]/N <= 0.1
 3 neff_ratio mu        0.621 high   N[eff]/N <= 0.1
 4 neff_ratio theta[1]  0.775 high   N[eff]/N <= 0.1
 5 neff_ratio eta[2]    0.777 high   N[eff]/N <= 0.1
 6 neff_ratio eta[5]    0.859 high   N[eff]/N <= 0.1
 7 neff_ratio eta[3]    0.897 high   N[eff]/N <= 0.1
 8 neff_ratio eta[1]    0.912 high   N[eff]/N <= 0.1
 9 neff_ratio theta[3]  0.915 high   N[eff]/N <= 0.1
10 neff_ratio theta[8]  0.917 high   N[eff]/N <= 0.1
11 neff_ratio eta[8]    0.919 high   N[eff]/N <= 0.1
12 neff_ratio eta[7]    0.979 high   N[eff]/N <= 0.1
13 neff_ratio eta[4]    0.984 high   N[eff]/N <= 0.1
14 neff_ratio eta[6]    0.991 high   N[eff]/N <= 0.1
15 neff_ratio theta[6]  0.998 high   N[eff]/N <= 0.1
16 neff_ratio theta[7]  1.02  high   N[eff]/N <= 0.1
17 neff_ratio theta[5]  1.02  high   N[eff]/N <= 0.1
18 neff_ratio theta[2]  1.03  high   N[eff]/N <= 0.1
19 neff_ratio theta[4]  1.03  high   N[eff]/N <= 0.1

mcmc_neff_hist(neff_ratio(fit))  ## stan_ess()形式で出力

ポイントは、neff_ratio()を通してから、出力させること。
~~stan_ess()でええやん~~

事後分布を確認する

ここからは、事後分布を様々な視点から確認していきたいと思います。

## デフォルト
stan_dens(fit, pars = c("mu", "tau"))

stan_dens(fit, pars = c("mu", "tau"), separate_chains = T)

# stan_dens(fit,pars=c('mu','tau'))$data

## beyasplot
mcmc_dens(samples, pars = c("mu", "tau"))

## chainごとに出力(ただ、一度に出す)
mcmc_dens_chains(samples, pars = c("mu", "tau"))

mcmc_dens_chains_data(samples, pars = c("mu", "tau"))  ## 上図のdataを取ってきたい場合

# A tibble: 8,192 x 6
   chain parameter interval_width     x density scaled_density
   <fct> <fct>              <dbl> <dbl>   <dbl>          <dbl>
 1 1     mu                     1 -6.52 0.00146         0.0174
 2 1     mu                     1 -6.49 0.00149         0.0178
 3 1     mu                     1 -6.46 0.00152         0.0182
 4 1     mu                     1 -6.42 0.00155         0.0185
 5 1     mu                     1 -6.39 0.00159         0.0189
 6 1     mu                     1 -6.36 0.00162         0.0193
 7 1     mu                     1 -6.33 0.00165         0.0197
 8 1     mu                     1 -6.30 0.00168         0.0201
 9 1     mu                     1 -6.27 0.00172         0.0205
10 1     mu                     1 -6.23 0.00175         0.0209
# ... with 8,182 more rows

## stan_dens(separate_chains=T)と同じようなこと実行
mcmc_dens_overlay(samples, pars = c("mu", "tau"))

ポイント、chian毎に出力するなら2パターンあって、片方は超絶醜い.
~~stan_densでええやん~~

複数のパラメータを一気に見たい場合

## デフォルト
stan_plot(fit, pars = "eta", ci_level = 0.95, outer_level = 1, show_density = T, 
    point_est = "mean")

## beyasplot
mcmc_areas(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, point_est = "mean")

mcmc_areas_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")

# A tibble: 16,524 x 5
   parameter interval interval_width     x density
   <fct>     <chr>             <dbl> <dbl>   <dbl>
 1 eta[1]    inner              0.95 -1.54  0.0575
 2 eta[1]    inner              0.95 -1.53  0.0578
 3 eta[1]    inner              0.95 -1.53  0.0581
 4 eta[1]    inner              0.95 -1.53  0.0584
 5 eta[1]    inner              0.95 -1.52  0.0588
 6 eta[1]    inner              0.95 -1.52  0.0591
 7 eta[1]    inner              0.95 -1.52  0.0594
 8 eta[1]    inner              0.95 -1.51  0.0598
 9 eta[1]    inner              0.95 -1.51  0.0601
10 eta[1]    inner              0.95 -1.51  0.0604
# ... with 16,514 more rows

## それぞれの事前分布が重なって表示されるようにする
## 図がコンパクトになるのが利点(公式)
mcmc_areas_ridges(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1)

mcmc_areas_ridges_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1)

# A tibble: 16,384 x 5
   parameter interval interval_width     x density
   <fct>     <chr>             <dbl> <dbl>   <dbl>
 1 eta[1]    inner              0.95 -1.54  0.0575
 2 eta[1]    inner              0.95 -1.53  0.0578
 3 eta[1]    inner              0.95 -1.53  0.0581
 4 eta[1]    inner              0.95 -1.53  0.0584
 5 eta[1]    inner              0.95 -1.52  0.0588
 6 eta[1]    inner              0.95 -1.52  0.0591
 7 eta[1]    inner              0.95 -1.52  0.0594
 8 eta[1]    inner              0.95 -1.51  0.0598
 9 eta[1]    inner              0.95 -1.51  0.0601
10 eta[1]    inner              0.95 -1.51  0.0604
# ... with 16,374 more rows

## 区間だけ図示したい
mcmc_intervals(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")

mcmc_intervals_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")

# A tibble: 8 x 9
  parameter outer_width inner_width point_est    ll     l        m     h
  <fct>           <dbl>       <dbl> <chr>     <dbl> <dbl>    <dbl> <dbl>
1 eta[1]              1        0.95 mean      -3.23 -1.54  0.399    2.18
2 eta[2]              1        0.95 mean      -3.32 -1.86 -0.00396  1.90
3 eta[3]              1        0.95 mean      -3.46 -2.14 -0.200    1.74
4 eta[4]              1        0.95 mean      -2.88 -1.82 -0.0575   1.65
5 eta[5]              1        0.95 mean      -3.51 -2.14 -0.361    1.50
6 eta[6]              1        0.95 mean      -3.38 -1.98 -0.227    1.65
7 eta[7]              1        0.95 mean      -3.13 -1.45  0.325    1.97
8 eta[8]              1        0.95 mean      -2.78 -1.78  0.0743   1.94
# ... with 1 more variable: hh <dbl>

ポイントは、regax_parsで正規表現が指定できるようになります。
今回のように“eta”だけ指定したい場合、“theta”もとってきてしまうので、“eta[”と指定して“eta”のみ抽出できるようにしているのが“^eta\[”です。

これに関しては、正規表現も簡単に使えるのが結構便利かも。

そして、最後は、事後分布(事後乱数)をヒストグラムやバイオリンプロット、散布図で書く場合

## デフォルト
stan_hist(fit, pars = c("mu", "tau"))

## beyasplot
mcmc_hist(samples, pars = c("mu", "tau"))

mcmc_hist_by_chain(samples, pars = c("mu", "tau"))  ## chain毎

## violin plot
mcmc_violin(samples, pars = "mu")

## 散布図
mcmc_scatter(samples, pars = c("eta[1]", "eta[2]"))

## ヒストグラム(density)と散布図同時に
mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("hist"))  ## 'hist' → 'dens'で変更可能

これは割とそのままですね～
mcmc_pairsは3変数以上でも出せますが、重くなるので注意が必要です。

NUTS診断系

以下は、HMCの細かい設定に関する出力や対数事後確率を出力したいときの関数になります。
出力が長いものは、head()で出力数を制限しています。

## デフォルト
stan_diag(fit)

## beyasplot

## nuts関連の値
nuts_params(fit) %>% head()

  Iteration     Parameter     Value Chain
1         1 accept_stat__ 0.9915962     1
2         2 accept_stat__ 0.9542163     1
3         3 accept_stat__ 0.9619325     1
4         4 accept_stat__ 0.9974984     1
5         5 accept_stat__ 0.9948058     1
6         6 accept_stat__ 0.6023858     1

## 対数事後確率
log_posterior(fit) %>% head()

  Iteration     Value Chain
1         1 -41.82965     1
2         2 -36.26125     1
3         3 -40.15663     1
4         4 -38.28819     1
5         5 -37.64542     1
6         6 -40.68307     1

## 上記を図で出力
mcmc_nuts_acceptance(nuts_params(fit), log_posterior(fit))

mcmc_nuts_divergence(nuts_params(fit), log_posterior(fit))

mcmc_nuts_energy(nuts_params(fit))

mcmc_nuts_stepsize(nuts_params(fit), log_posterior(fit))

mcmc_nuts_treedepth(nuts_params(fit), log_posterior(fit))

並行でだす関数(mcmc_parcoord)

複数のパラメータを同時に(並行に)出力する関数です。
条件間・参加者間を通して確認したいときに便利かもしれません。
ただ、出力が結構重いので注意が必要です。

## bayesplot
mcmc_parcoord(samples, regex_pars = "^eta")

mcmc_parcoord_data(samples, regex_pars = "^eta") %>% head

# A tibble: 6 x 4
   Draw Parameter    Value Divergent
  <int> <fct>        <dbl>     <dbl>
1     1 eta[1]    -0.0968          0
2     2 eta[1]     0.736           0
3     3 eta[1]     1.25            0
4     4 eta[1]     1.30            0
5     5 eta[1]    -0.00446         0
6     6 eta[1]    -0.202           0

パラメータリカバリ

モデルの確認のためにパラメータリカバリの性能を確かめることは多々あるとは思いますが、真値と推定値の確認に便利な関数が以下になります。

## 真値データの作成(8shcool.stan用)

J <- 10

set.seed(1234)
eta <- rnorm(J, 0, 1)
mu <- 7.5  ## 適当です
tau <- 6.3  ## 適当です 
sigma <- abs(rnorm(J, 0, 1))
theta <- NULL
y <- NULL
for (j in 1:J) {
    theta[j] <- mu + tau * eta[j]
    y[j] <- rnorm(1, theta[j], sigma[j])
}
true <- c(mu, tau, eta, theta)

## Stan data 指定
schools_dat <- list(J = J, y = y, sigma = sigma)

## stan を回す(Modelを読み込んでいる前提です。)
fit <- sampling(object = model, data = schools_dat, iter = 1000, chains = 4)

## 結果 ##
fit

Inference for Stan model: 4a38cae9bad394789360bce03f17c3ed.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff
mu          4.90    0.17 2.71  -0.39   3.30   4.84   6.43  10.13   250
tau         7.46    0.14 2.22   4.53   5.86   6.96   8.53  13.32   244
eta[1]     -0.70    0.02 0.38  -1.43  -0.95  -0.71  -0.45   0.06   329
eta[2]      0.55    0.02 0.40  -0.21   0.29   0.55   0.82   1.35   306
eta[3]      1.30    0.03 0.51   0.34   0.96   1.29   1.62   2.33   272
eta[4]     -1.75    0.03 0.55  -2.84  -2.12  -1.75  -1.38  -0.65   301
eta[5]      0.66    0.02 0.41  -0.13   0.39   0.65   0.93   1.49   299
eta[6]      0.82    0.03 0.42   0.04   0.54   0.82   1.10   1.67   271
eta[7]     -0.10    0.02 0.35  -0.78  -0.33  -0.10   0.14   0.61   313
eta[8]     -0.25    0.02 0.36  -0.97  -0.48  -0.25  -0.01   0.46   359
eta[9]     -0.13    0.02 0.35  -0.82  -0.37  -0.13   0.09   0.57   323
eta[10]    -0.65    0.02 0.46  -1.58  -0.96  -0.65  -0.34   0.25   487
theta[1]   -0.01    0.01 0.49  -0.94  -0.32  -0.01   0.32   0.95  2289
theta[2]    8.68    0.02 1.00   6.71   8.01   8.69   9.34  10.59  2247
theta[3]   13.87    0.01 0.79  12.33  13.34  13.87  14.39  15.45  3180
theta[4]   -7.25    0.00 0.07  -7.37  -7.29  -7.25  -7.20  -7.12  1910
theta[5]    9.44    0.02 0.93   7.62   8.80   9.44  10.09  11.25  2310
theta[6]   10.53    0.00 0.11  10.31  10.45  10.52  10.60  10.76  1787
theta[7]    4.17    0.01 0.52   3.19   3.83   4.16   4.52   5.20  2159
theta[8]    3.14    0.02 0.91   1.37   2.52   3.13   3.74   4.86  2065
theta[9]    3.94    0.02 0.82   2.27   3.39   3.92   4.46   5.54  2074
theta[10]   0.30    0.07 2.32  -4.02  -1.29   0.27   1.83   5.28   968
lp__      -19.94    0.18 3.30 -26.80 -22.05 -19.72 -17.41 -14.41   328
          Rhat
mu        1.01
tau       1.01
eta[1]    1.01
eta[2]    1.01
eta[3]    1.01
eta[4]    1.01
eta[5]    1.01
eta[6]    1.01
eta[7]    1.01
eta[8]    1.01
eta[9]    1.01
eta[10]   1.01
theta[1]  1.00
theta[2]  1.00
theta[3]  1.00
theta[4]  1.00
theta[5]  1.00
theta[6]  1.00
theta[7]  1.00
theta[8]  1.00
theta[9]  1.00
theta[10] 1.00
lp__      1.00

Samples were drawn using NUTS(diag_e) at Sat Dec 01 12:43:12 2018.
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).

## true(真値)との比較
matfit <- as.matrix(fit)

## 推定値と真値の比較
mcmc_recover_intervals(matfit[, -ncol(matfit)], true)  ## fitの最後にでてくる'lp__'の分を抜いている

## ヒストグラム(vs 真値)
mcmc_recover_hist(matfit[, 2:3], true[2:3])

## 散布図(vs 真値) だとおもうのですがこれだけうまくいきません。

# mcmc_recover_scatter(matfit[,-ncol(matfit)], true,point_est = 'mean')

真値と簡単に比較できるのは、便利ですね。
ただ、trueの入れる順序を間違うと誤った解釈につながる可能性ありますし、そこら辺注意ですね・・・。

関数を探したいときに便利な関数

これだけ多くのmcmc_関数があるととても便利な反面、いざ使うと気になると思い出せないですよね。
かといって、暗記するのも大変だし…

そんなときに便利なのが、関数を探すための関数があります。

## bayesplot
available_mcmc()  ## mcmc_関数の種類を提示

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_data
  mcmc_areas_ridges
  mcmc_areas_ridges_data
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_chains_data
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_intervals_data
  mcmc_neff
  mcmc_neff_data
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_parcoord_data
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_data
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

## 特定の用途に関する関数を探したい(ここでは、rhatもしくはn_eff関連の関数を検索)。
available_mcmc("rhat|neff")

bayesplot MCMC module:
(matching pattern 'rhat|neff') 
  mcmc_neff
  mcmc_neff_data
  mcmc_neff_hist
  mcmc_rhat
  mcmc_rhat_data
  mcmc_rhat_hist

とりあえずこのパッケージで遊びたい

そんな方のために、デモデータおよびデモ事後乱数も付いています。
こちらも参考まで

## example_ シリーズ

example_mcmc_draws() %>% head()

[1] -14.147624 -20.037724 -20.982792 -36.288584  -7.581124 -10.409814

example_group_data() %>% head()

[1] GroupB GroupB GroupB GroupB GroupB GroupA
Levels: GroupA GroupB

example_x_data() %>% head()

[1] 121.11753  89.36188 115.44316  99.44964  92.74571 107.90184

example_y_data() %>% head()

[1]  65  98  85  83 115  98

# example_yrep_draws() %>% head() ##
# 表示数が多いのでコメントアウトしています。

ここらへんを使えば、ppc_シリーズの挙動も簡単に調べられると思います。
ただし、ppc_シリーズは、少しクセが強かったり、パラメータ数増やすと結構見にくかったりで大変なので頑張ってくださいｗ

まとめ

こんなかんじでbayesplotの特にmcmc_シリーズを見てきました。
便利な部分と既存の関数が良い部分と様々でした。

自分の用途や余裕に合わせて使いこなせたらいいですね。Enjoy!!

今回使用した全Rコード

# setwd('')

rm(list = ls())
library(tidyverse)  ## ほぼパイプのために読み込んでいます。
library(rstan)
library(bayesplot)  ## こいつが今回の主役です。
rstan_options(auto_write = TRUE)  ## なくても構いません
options(mc.cores = parallel::detectCores())  ## なくても構いません
options(max.print = 99999)  ## なくても構いません

schools_dat <- list(J = 8, y = c(28, 8, -3, 7, -1, 1, 18, 12), sigma = c(15, 
    10, 16, 11, 9, 11, 10, 18))

model <- stan_model("model.stan")

fit <- sampling(object = model, data = schools_dat, iter = 1000, chains = 4)

## 結果 ##
fit

## デフォルト not bayesplot
stan_rhat(fit)
stan_rhat(fit)$data  ## dataとしてとってきたい

## bayesplot

## rhatがdataとして出力される
rhat(fit)

## rhatが図として出力
mcmc_rhat(rhat(fit))

## 表示した図に関するdata
mcmc_rhat_data(rhat(fit))

## histgramとして表示(stan_rhatに近い形で出力)
mcmc_rhat_hist(rhat(fit))

## 特定のパラメータのR hat をみたい。
stan_rhat(fit, pars = c("mu", "eta"))  ## stanの関数
mcmc_rhat(rhat(fit, pars = c("mu", "eta")))  ## bayesplotの関数

## デフォルト
stan_trace(fit, inc_warmup = T)  ## inc_warmupでwarmup区間指定可能
traceplot(fit, inc_warmup = T)

## bayesplot
samples <- as.array(fit)
mcmc_trace(samples, pars = "mu")  ## n_warmupでwarmup区間指定可能

## 特定のchainを特に見たい(目立たせて確認したい)
mcmc_trace_highlight(samples, pars = "mu", highlight = 2)

## 事後分布と同時に確認したい(パラメータ多いと重いので指定推奨)
mcmc_combo(samples, pars = c("mu"))

## デフォルト
stan_ac(fit, pars = c("mu", "tau"))

## chain別
stan_ac(fit, pars = c("mu", "tau"), separate_chains = T)

## bayesplot

## 折れ線表記
mcmc_acf(samples, pars = c("mu", "tau"))

## 棒グラフ表記
mcmc_acf_bar(samples, pars = c("mu", "tau"))

## デフォルト
stan_ess(fit)
stan_ess(fit)$data  ## dataとして出力

## bayesplot n_effの値をdataとして確認
neff_ratio(fit)

## 図で確認
mcmc_neff(neff_ratio(fit))
mcmc_neff_data(neff_ratio(fit))  ## 図で出力したdataを出力
mcmc_neff_hist(neff_ratio(fit))  ## stan_ess()形式で出力

## デフォルト
stan_dens(fit, pars = c("mu", "tau"))
stan_dens(fit, pars = c("mu", "tau"), separate_chains = T)
# stan_dens(fit,pars=c('mu','tau'))$data

## beyasplot
mcmc_dens(samples, pars = c("mu", "tau"))

## chainごとに出力(ただ、一度に出す)
mcmc_dens_chains(samples, pars = c("mu", "tau"))
mcmc_dens_chains_data(samples, pars = c("mu", "tau"))  ## 上図のdataを取ってきたい場合

## stan_dens(separate_chains=T)と同じようなこと実行
mcmc_dens_overlay(samples, pars = c("mu", "tau"))

## デフォルト
stan_plot(fit, pars = "eta", ci_level = 0.95, outer_level = 1, show_density = T, 
    point_est = "mean")

## beyasplot
mcmc_areas(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, point_est = "mean")
mcmc_areas_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")

## それぞれの事前分布が重なって表示されるようにする
## 図がコンパクトになるのが利点(公式)
mcmc_areas_ridges(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1)
mcmc_areas_ridges_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1)

## 区間だけ図示したい
mcmc_intervals(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")
mcmc_intervals_data(samples, regex_pars = "^eta\\[", prob = 0.95, prob_outer = 1, 
    point_est = "mean")

## デフォルト
stan_hist(fit, pars = c("mu", "tau"))

## beyasplot
mcmc_hist(samples, pars = c("mu", "tau"))
mcmc_hist_by_chain(samples, pars = c("mu", "tau"))  ## chain毎

## violin plot
mcmc_violin(samples, pars = "mu")

## 散布図
mcmc_scatter(samples, pars = c("eta[1]", "eta[2]"))

## ヒストグラム(density)と散布図同時に
mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("hist"))  ## 'hist' → 'dens'で変更可能

## デフォルト
stan_diag(fit)

## beyasplot

## nuts関連の値
nuts_params(fit) %>% head()

## 対数事後確率
log_posterior(fit) %>% head()

## 上記を図で出力
mcmc_nuts_acceptance(nuts_params(fit), log_posterior(fit))
mcmc_nuts_divergence(nuts_params(fit), log_posterior(fit))
mcmc_nuts_energy(nuts_params(fit))
mcmc_nuts_stepsize(nuts_params(fit), log_posterior(fit))
mcmc_nuts_treedepth(nuts_params(fit), log_posterior(fit))

## bayesplot
mcmc_parcoord(samples, regex_pars = "^eta")
mcmc_parcoord_data(samples, regex_pars = "^eta") %>% head

## 真値データの作成(8shcool.stan用)

J <- 10

set.seed(1234)
eta <- rnorm(J, 0, 1)
mu <- 7.5  ## 適当です
tau <- 6.3  ## 適当です 
sigma <- abs(rnorm(J, 0, 1))
theta <- NULL
y <- NULL
for (j in 1:J) {
    theta[j] <- mu + tau * eta[j]
    y[j] <- rnorm(1, theta[j], sigma[j])
}
true <- c(mu, tau, eta, theta)

## Stan data 指定
schools_dat <- list(J = J, y = y, sigma = sigma)

## stan を回す(Modelを読み込んでいる前提です。)
fit <- sampling(object = model, data = schools_dat, iter = 1000, chains = 4)

## 結果 ##
fit

## true(真値)との比較
matfit <- as.matrix(fit)

## 推定値と真値の比較
mcmc_recover_intervals(matfit[, -ncol(matfit)], true)  ## fitの最後にでてくる'lp__'の分を抜いている

## ヒストグラム(vs 真値)
mcmc_recover_hist(matfit[, 2:3], true[2:3])

## 散布図(vs 真値) だとおもうのですがこれだけうまくいきません。

# mcmc_recover_scatter(matfit[,-ncol(matfit)], true,point_est = 'mean')

## bayesplot
available_mcmc()  ## mcmc_関数の種類を提示

## 特定の用途に関する関数を探したい(ここでは、rhatもしくはn_eff関連の関数を検索)。
available_mcmc("rhat|neff")

## example_ シリーズ
example_mcmc_draws() %>% head()
example_group_data() %>% head()
example_x_data() %>% head()
example_y_data() %>% head()
example_yrep_draws() %>% head()