StanのMCMC結果&パラメータ結果を可視化する① bayesplotとShinyStan
背景
最近rstan経由でStanを使ってる。rstanを用いた結果(収束診断とか事後確率分布とか)はそのままのデータでは可視化をするのが面倒。
可視化するのに便利なパッケージはないか調べてみると、ggmcmc
とか bayesplot
とか shinystan
とか tidybayes
とか色々ある模様。
友人に使い分けを聞いたりした感じだと、
- MCMCがうまくいってるかとか見るなら
bayesplot
かshinystan
- パラメータの事後分布(サンプリング結果)は
tidybayes
という使い分けがいいそうな。 とりあえずそれぞれで試してみる。
準備
データは以下のサイトと同様のものを使用する。
Introduction to bayesplot (mcmc_ series)
bayesplotでMCMC結果を見る
内容も上記サイトをトレースする。
概要
bayesplotではrstanの標準関数より便利な形でggplotオブジェクトとして結果を返してくれる。
bayesplotパッケージでは、大きく二つの関数系に分かれているとのこと - 1. mcmc xxx(収束や事後分布を確認するための関数群) - 2. ppc xxx(事後予測チェックを行うのに便利な関数群)
パラメータ推定が収束しているか
の確認
rhat(fit)
で書くパラメータの 結果を出力できる。
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を出力してくれる。
は一般的に1.1以下なら良いため、その基準に基づいて凡例が色分けされている。
mcmc_rhat_hist(rhat(fit))
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"))
パラメータ推定が安定しているか
自己相関の確認
先ほどと同様に、配列型を関数に渡す。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
実行有効サンプルサイズ(ESS(n_eff))
0.1以下だと自己相関が高く、MCMCが上手くいってないとのこと(あんまよく理解してない
neff_ratio()
でn_eff値を数値で、mcmc_neff()
でn_eff値を図として確認できる。
# 値 neff_ratio(fit) # 図で出力 mcmc_neff_hist(neff_ratio(fit))
事後分布の確認
冒頭で書いたように、tidybayes
の方が良さげだが一応。
確率密度
# chain合算 g3 = mcmc_dens(samples, pars = c("mu", "tau")) # chain毎 g4 = mcmc_dens_overlay(samples, pars = c("mu", "tau")) g3 / g4
まとめて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")
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"))
chain毎
mcmc_hist_by_chain(samples, pars = c("mu", "tau"))
散布図
mcmc_scatter(samples, pars = c("mu", "tau"))
散布図とヒストグラムまとめて
mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("hist")) mcmc_pairs(samples, pars = c("mu", "tau"), diag_fun = c("dens"))
ShinyStanでMCMC結果を見る
ShinyStanは名前通りShinyで動いている。
要はbayesplotとかで出力していることを、まとめてShinyアプリ上で勝手にplotしてくれるという優れもの。
rstanの結果のほかに、rstanarm
やbrms
のような、ベイズ系パッケージに対しても使えるらしい。
グラフ画像を使ってなにかしたり、より自分好みの画像を作るとかじゃない限りは基本的にShinyStanを使えばいい気がする。
以下のように、作成したモデルに`launch_shinystan()'関数を適用すると起動する。
library(shinystan) my_sso = launch_shinystan(fit)
メニューはそれぞれ
- DIAGNOSE:MCMCの収束診断
- ESTIMATE:パラメータの推定結果
- EXPLORE:グラフ
- MORE:その他(コードや単語解説など)
となっている。
少し表示内容が変わっているが内容の詳細は以下のP37あたりを参照するとよい。
tidybayes
長くなったのでtidybayesは次の記事で。
参考
Plotting for Bayesian Models • bayesplot