Bayesplotパッケージはいいよ

ベイズの話があちこちで盛り上がっていますね!

ところで,ベイジアン・モデリングで得られた結果を表現する切り口っていろいろあるのですが,これを一手に引き受けてくれる便利なパッケージ,bayespotの存在を知りました。

これまではshinystanでみると綺麗!だったのですが,ブラウザが立ち上がるのでちょっと面倒。bayesplotパッケージはplot領域に出力されるし,結果がggplot2オブジェクトなので,ggplot2のレイアを上書きして言ったりすることができるのも便利。

ということで,練習がてら挙動を確認してみる。

まずは簡単なモデルを書いてみることに。

# ライブラリの読み込み
library(rstan)
##  要求されたパッケージ ggplot2 をロード中です
##  要求されたパッケージ StanHeaders をロード中です
## rstan (Version 2.14.1, packaged: 2016-12-28 14:55:41 UTC, GitRev: 5fa1e80eb817)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(bayesplot)
## This is bayesplot version 1.1.0
# 擬似データ作成
N <- 100
x <- runif(N,1,10)
yhat <- 7 + 15 * x 
y <- yhat + rnorm(N,0,10)


# Stanコード
stancode <- '
data{
  int<lower=0> N;
  vector[N] X;
  vector[N] Y;
}

parameters{
  real mu;
  real beta;
  real<lower=0> sigma;
}

model{
  Y ~ normal(mu+beta*X,sigma);
}

generated quantities{
  vector[N] pred;
  for(n in 1:N)
    pred[n] = normal_rng(mu+beta*X[n],sigma);
}
'

# モデルフィット
fit <- stan(model_code = stancode,data=list(N=N,X=x,Y=y))
print(fit,pars=c("mu","beta","sigma"))
## Inference for Stan model: de83821dfd3dddca7f339079f32d2db5.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## mu     4.83    0.06 2.14  0.69  3.41  4.86  6.25  9.08  1488    1
## beta  15.29    0.01 0.34 14.61 15.07 15.28 15.51 15.97  1489    1
## sigma  9.76    0.01 0.70  8.50  9.28  9.71 10.20 11.20  2317    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Mar 24 15:00:38 2017.
## 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パッケージに与えて色々するんだけど,stanfitオブジェクトのままではうまくいかないので,MCMCサンプルをarray 型にしておく。これでもって,サンプリングのいろいろな結果を描いてみる。

fit.array <- as.array(fit)

# MCMCのトレースプロット
mcmc_trace(fit.array,pars=c("mu","beta","sigma"))

plot of chunk unnamed-chunk-2

# 自己相関
mcmc_acf(fit.array,pars=c("mu","beta","sigma"))

plot of chunk unnamed-chunk-2

# 事後確率密度
mcmc_dens(fit.array,pars=c("mu","beta","sigma"))

plot of chunk unnamed-chunk-2

mcmc_dens_overlay(fit.array,pars=c("mu","beta","sigma"))

plot of chunk unnamed-chunk-2

# ヒストグラム
mcmc_hist(fit.array,pars=c("mu","beta","sigma"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot of chunk unnamed-chunk-2

mcmc_hist_by_chain(fit.array,pars=c("mu","beta","sigma"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot of chunk unnamed-chunk-2

# ヴァイオリンプロット
mcmc_violin(fit.array,pars=c("mu","beta","sigma"))

plot of chunk unnamed-chunk-2

# 区間推定プロット
mcmc_areas(fit.array,pars=c("mu","beta","sigma"),prob=0.95)

plot of chunk unnamed-chunk-2

収束判断に使うRhatやn_effは関数で呼び出して図示する。

mcmc_rhat(rhat(fit))

plot of chunk conversion

mcmc_neff(neff_ratio(fit))

plot of chunk conversion

いいな,と思ったのは事後予測チェックのプロットもできること。generated quantitiesで推定値から事後予測用のデータを発生させて置いたのはこのためです。

# 生成量としての予測値を抜き出す
Yrep <- rstan::extract(fit,pars="pred")$pred
# 事後予測チェックPosterior Predictive Checking の一連の関数
# 
ppc_stat(y,Yrep,stat="mean")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot of chunk unnamed-chunk-3

ppc_stat_2d(y,Yrep,stat=c("mean","sd"))

plot of chunk unnamed-chunk-3

いくつかのサンプル(ここでは最初の9つ)とデータを並べてみる。
サンプルを指定しないと(4000あるので)反応が返ってこなくなっちゃう!(ヘルプにも「小さい数字でね」って書いてました)

ppc_dens(y,Yrep[1:9,])

plot of chunk unnamed-chunk-4

ppc_dens_overlay(y,Yrep[1:9,])

plot of chunk unnamed-chunk-4

ppc_boxplot(y,Yrep[1:9,])

plot of chunk unnamed-chunk-4

ppc_freqpoly(y,Yrep[1:9,])
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.