NagoyaStat#8 アヒル本5章
あっ、どーも僕です。
NagoyaStat
昨日、NagoyaStat#8が開催され、スピーカーとして参加しました。 今回のお題はアヒル本こと『StanとRでベイズ統計モデリング』の5章です。 ここでは発表資料は公開します。よろしければ参照ください。
なお、資料は本ありきの構成です。 インターネットで探すと、自分のデータに適用したりと、 本以上のことをやられているかたもおられましたが、 今回はそこまでできませんでした。
ユーティリティの作成
5章の目的の1つはいまのうちに必要な関数をつくることだと感じました。 そこで、ここでは私が作ったものを一部載せます。
Stanのコンパイル
私だけかも知れませんが、rmdを書きながらStanを実行すると やたらとRが落ちました。 そこで安全1にStanのコンパイルとサンプリングができるラッパーを用意しました。
作成した関数は2つです。create_samplerが安全に実行するための関数です。
- create_sampler:モデルファイルのパスを引数にして、sampling関数を部分評価したものを返す
- summary2:summary(stanfit)$summaryをして、データフレームで返す
例えばアヒル本5章の二項ロジスティック回帰を行う際には、 私は下記のスクリプトのようにしています。
data <- attend2 %>% as.list() %>% update_list( PersonID = NULL, Score = ~ Score/ 200, N = ~ length( Y ) ) model_file <- "src/model-5-4.stan" compiled_sampler <- create_sampler(model_file, save = TRUE) fit <- compiled_sampler( data = data, seed = 1234, chain = 1 ) smr <- summary2(fit); write_csv( smr, paste0(model_file, ".summary.csv" ) )
create_samplerの中身は次のようになっています。 やっていることとしてはコンパイル済みのStanオブジェクトがあればそれを読み込み、 無ければコンパイルしてから読み込むというものです。
ポイントはcompile_stan_file の最後の1行です。 この1行ではコンパイル済みのStanオブジェクトを持つRオブジェクトを削除して、 ガーベージコレクションしています。 わたしの場合、コンパイルをしばかりのホヤホヤのオブジェクトを持っていると、 Rが不安定になっていました2
compile_stan_file <- function (stan_file) { ofile <- paste0( stan_file, ".rds" ) cmodel <- stan_model( file = stan_file ) write_rds( cmodel, ofile ) rm( cmodel ); gc(); gc() } read_compiled_stan_file <- function (stan_file, compile = FALSE) { compiled_stan_file <- paste0( stan_file, ".rds" ) if (!file.exists( compiled_stan_file ) || compile ) compile_stan_file( stan_file ) return( read_rds( compiled_stan_file ) ) } generate_opath <- function (dir, file, ext) { file.path( dir, paste0( file, ext ) ) } create_sampler <- function (stan_file, compile = FALSE, save = TRUE, save_dir = dirname( stan_file ) ) { args <- list( ...f = sampling, .lazy = FALSE, object = read_compiled_stan_file( stan_file, compile ), sample_file = generate_opath( save_dir, basename( stan_file ), ".sample.csv" ), diagnostic_file = generate_opath( save_dir , basename( stan_file ), ".diagnostic.csv" ) ) if (!save) { args <- list_modify( args, sample_file = NULL, diagnostic_file = NULL ) } if (!dir.exists( save_dir )) { dir.create( save_dir ) } spr <- lift( partial )( args ) spr }
summary2の中身は次です。
mat2tbl <- function (mat, rname = NULL) { as_data_frame( mat ) %>% mutate( varname = rname ) %>% select( varname, everything() ) } summary2 <- compose( partial( mat2tbl, rname = rownames( ... ) ), partial( pluck, "summary", .first = FALSE ), summary )
予測値の分布
5章では回帰係数にせよ、応答変数にせよ、MCMCサンプルの10%、50%、90%タイル値を使って グラフを描くことが多かったです。そこで、MCMCサンプルを引数に タイル値を計算し、列に追加してくれる関数をつくりました。 これについてはもっと良さげな関数が提供されていると思いましたが、 調べた限りではわからなかったです。
で、その関数というのがこちらです。
q10 <- partial( quantile, probs = .1 ) q50 <- partial( quantile, probs = .5 ) q90 <- partial( quantile, probs = .9 ) qqs <- function( ms ) { partial( mutate, q10 = map_dbl(ms, q10), q50 = map_dbl(ms, q50), q90 = map_dbl(ms, q90), .first = FALSE ) }
例えば、5.2の二項ロジスティック回帰で回帰係数の分布を描く際には、次のような感じで使っています。 抽出したサンプリング結果を上記の関数に与えてできる関数に、 データフレームを与えることで、タイル値をそのデータフレームにくっつけてくれます。 後はggplotで簡単に作図できます。
pars <- c("b1", "b2", "b3") ms <- extract( fit, pars = pars) ret <- data_frame( x = seq.int(length(pars)) ) %>% qqs(ms)() ret %>% ggplot() + geom_linerange(aes(x, ymin = q10, ymax = q90)) + geom_point(aes(x, q50), size = 5) + theme_bw()
お礼
会場は今回もヤフー株式会社様から提供いただきました。 いつもありがとうございます。 また、無理な日程を承知してくださった NagoyaStatの皆様ありがとうございました。
誰かと自分の未来のために