まずは蝋の翼から。

学んだことを書きながら確認・整理するためのメモブログ。こういうことなのかな?といったことをふわっと書いたりしていますが、理解が浅いゆえに的はずれなことも多々あると思うのでツッコミ歓迎

StanのMCMC結果&パラメータ結果を可視化する① bayesplotとShinyStan

背景

最近rstan経由でStanを使ってる。rstanを用いた結果(収束診断とか事後確率分布とか)はそのままのデータでは可視化をするのが面倒。
可視化するのに便利なパッケージはないか調べてみると、ggmcmc とか bayesplot とか shinystan とか tidybayes とか色々ある模様。

友人に使い分けを聞いたりした感じだと、

  • MCMCがうまくいってるかとか見るならbayesplotshinystan
  • パラメータの事後分布(サンプリング結果)はtidybayes

という使い分けがいいそうな。 とりあえずそれぞれで試してみる。

準備

データは以下のサイトと同様のものを使用する。

Introduction to bayesplot (mcmc_ series)

bayesplotでMCMC結果を見る

内容も上記サイトをトレースする。

概要

bayesplotではrstanの標準関数より便利な形でggplotオブジェクトとして結果を返してくれる。

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

パラメータ推定が収束しているか

\hat{R} の確認

rhat(fit)で書くパラメータの \hat{R} 結果を出力できる。

rhat(fit)

# =>
       mu       tau    eta[1]    eta[2]    eta[3]    eta[4]    eta[5]    eta[6]    eta[7] 
1.0112444 1.0086249 1.0002212 1.0000860 1.0018156 0.9996789 1.0001951 0.9985509 1.0004175 
   eta[8]  theta[1]  theta[2]  theta[3]  theta[4]  theta[5]  theta[6]  theta[7]  theta[8] 
1.0006213 1.0041166 0.9996360 1.0007257 1.0001193 1.0005308 1.0009002 1.0019892 1.0022765 
     lp__ 
1.0126542 

このデータをmcmc_rhat_hist()に渡すといい感じにrhatを出力してくれる。
\hat{R} は一般的に1.1以下なら良いため、その基準に基づいて凡例が色分けされている。

mcmc_rhat_hist(rhat(fit))

f:id:chito_ng:20190518205209p:plain

traceplotの確認

一度stanの結果をarrayにしてからmcmc_trace()にわたす。mcmc_combo()に渡した場合は事後分布と共に表示してくれる。

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

## 事後分布と同時(パラメータ多いと重い)
mcmc_combo(samples, pars = c("mu","tau"))

f:id:chito_ng:20190518205934p:plain

パラメータ推定が安定しているか

自己相関の確認

先ほどと同様に、配列型を関数に渡す。mcmc_acf()で折れ線、mcmc_acf_bar()で棒グラフとして、lagとの相関が(自己相関)が見れる。

library(patchwork)

## 折れ線表記
g1 = mcmc_acf(samples, pars = c("mu", "tau"))
## 棒グラフ表記
g2 = mcmc_acf_bar(samples, pars = c("mu", "tau"))

g1 | g2

f:id:chito_ng:20190518210315p:plain

実行有効サンプルサイズ(ESS(n_eff))

0.1以下だと自己相関が高く、MCMCが上手くいってないとのこと(あんまよく理解してない
neff_ratio()でn_eff値を数値で、mcmc_neff()でn_eff値を図として確認できる。

# 値
neff_ratio(fit)
# 図で出力
mcmc_neff_hist(neff_ratio(fit))

f:id:chito_ng:20190518211808p:plain

事後分布の確認

冒頭で書いたように、tidybayesの方が良さげだが一応。

確率密度

# chain合算
g3 = mcmc_dens(samples, pars = c("mu", "tau"))

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

g3 / g4

f:id:chito_ng:20190518212417p:plain

まとめて1グラフに表示したい場合。上記のmcmc_dens()でもひとつといえばひとつなので好みか?
ただ、こっちのmcmc_areas()の場合は正規表現が使えるようになったり、信頼区間出たりする。

mcmc_areas(samples, c("mu", "tau"), prob = 0.95, prob_outer = 1, point_est = "mean")

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

f:id:chito_ng:20190518212824p:plain

f:id:chito_ng:20190518212735p:plain

mcmc_areas_data()で値として出力できる。

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

# A tibble: 16,534 x 5
   parameter interval interval_width     x density
   <fct>     <chr>             <dbl> <dbl>   <dbl>
 1 eta[1]    inner              0.95 -1.47  0.0709
 2 eta[1]    inner              0.95 -1.46  0.0711
 3 eta[1]    inner              0.95 -1.46  0.0713
 4 eta[1]    inner              0.95 -1.45  0.0716
 5 eta[1]    inner              0.95 -1.45  0.0718
 6 eta[1]    inner              0.95 -1.45  0.0720
 7 eta[1]    inner              0.95 -1.44  0.0722
 8 eta[1]    inner              0.95 -1.44  0.0723
 9 eta[1]    inner              0.95 -1.44  0.0725
10 eta[1]    inner              0.95 -1.43  0.0727
# … with 16,524 more rows

ヒストグラム

chain合算

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

f:id:chito_ng:20190518220034p:plain

chain毎

mcmc_hist_by_chain(samples, pars = c("mu", "tau"))

f:id:chito_ng:20190518220153p:plain

散布図

mcmc_scatter(samples, pars = c("mu", "tau"))

f:id:chito_ng:20190518220359p:plain

散布図とヒストグラムまとめて

mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("hist"))

mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("dens"))

f:id:chito_ng:20190518220624p:plain

f:id:chito_ng:20190518220636p:plain

ShinyStanでMCMC結果を見る

ShinyStanは名前通りShinyで動いている。

要はbayesplotとかで出力していることを、まとめてShinyアプリ上で勝手にplotしてくれるという優れもの。 rstanの結果のほかに、rstanarmbrmsのような、ベイズ系パッケージに対しても使えるらしい。

Getting Started • shinystan

グラフ画像を使ってなにかしたり、より自分好みの画像を作るとかじゃない限りは基本的にShinyStanを使えばいい気がする。

以下のように、作成したモデルに`launch_shinystan()'関数を適用すると起動する。

library(shinystan)

my_sso = launch_shinystan(fit)

f:id:chito_ng:20190519002248p:plain

メニューはそれぞれ

  • DIAGNOSE:MCMCの収束診断
  • ESTIMATE:パラメータの推定結果
  • EXPLORE:グラフ
  • MORE:その他(コードや単語解説など)

となっている。
少し表示内容が変わっているが内容の詳細は以下のP37あたりを参照するとよい。

www.slideshare.net

tidybayes

長くなったのでtidybayesは次の記事で。

参考

Plotting for Bayesian Models • bayesplot

Introduction to bayesplot (mcmc_ series)

Bayesplotパッケージはいいよ – Kosugitti's BLOG