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"))

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

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

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

# ヒストグラム
mcmc_hist(fit.array,pars=c("mu","beta","sigma"))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

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

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

mcmc_rhat(rhat(fit))

mcmc_neff(neff_ratio(fit))

いいな，と思ったのは事後予測チェックのプロットもできること。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`.

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

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

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

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

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

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

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

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

二次元に並べる場合。軸の範囲はggplotのlims関数で指定している

ppc_scatter_avg(y,Yrep)+ggplot2::lims(x=c(min(y),max(y)),y=c(min(y),max(y)))

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

予測値とデータとの差＝errorのプロット

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

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

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

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

サンプルの平均を取ってみる場合

ppc_error_scatter_avg(y,Yrep)

独立変数Xと平均的な誤差

ppc_error_scatter_avg_vs_x(y,Yrep,x)

これらのプロットは自分で作ってもいいけど，こうやって関数化されていると楽だよねえ。
モデルがrstanarm経由で作られたりしているともっと簡単ですよ。

ということで今日はここまで。