ベイズの話があちこちで盛り上がっていますね!
ところで,ベイジアン・モデリングで得られた結果を表現する切り口っていろいろあるのですが,これを一手に引き受けてくれる便利なパッケージ,bayespotの存在を知りました。
これまではshinystanでみると綺麗!だったのですが,ブラウザが立ち上がるのでちょっと面倒。bayesplotパッケージはplot領域に出力されるし,結果がggplot2オブジェクトなので,ggplot2のレイアを上書きして言ったりすることができるのも便利。
ということで,練習がてら挙動を確認してみる。
まずは簡単なモデルを書いてみることに。
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);
}
# ライブラリの読み込み
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
library(bayesplot)
# 擬似データ作成
N <- 100
x <- runif(N,1,10)
yhat <- 7 + 15 * x
y <- yhat + rnorm(N,0,10)
# モデルフィット
fit <- sampling(sampleCode,data=list(N=N,X=x,Y=y))
print(fit,pars=c("mu","beta","sigma"))
## Inference for Stan model: 3c6ac53b227702d2d44ca516d9974187. ## 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 6.82 0.07 2.67 1.54 5.02 6.87 8.59 12.05 1506 1 ## beta 14.94 0.01 0.43 14.09 14.65 14.93 15.23 15.78 1510 1 ## sigma 11.67 0.02 0.86 10.15 11.06 11.62 12.22 13.46 2246 1 ## ## Samples were drawn using NUTS(diag_e) at Wed Nov 29 12:19:32 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,])