NagoyaStat

昨日、NagoyaStat#8が開催され、スピーカーとして参加しました。今回のお題はアヒル本こと『StanとRでベイズ統計モデリング』の5章です。ここでは発表資料は公開します。よろしければ参照ください。

なお、資料は本ありきの構成です。インターネットで探すと、自分のデータに適用したりと、本以上のことをやられているかたもおられましたが、今回はそこまでできませんでした。

ユーティリティの作成

5章の目的の1つはいまのうちに必要な関数をつくることだと感じました。そこで、ここでは私が作ったものを一部載せます。

Stanのコンパイル

私だけかも知れませんが、rmdを書きながらStanを実行するとやたらとRが落ちました。そこで安全¹に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が不安定になっていました²

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の皆様ありがとうございました。

誰かと自分の未来のために

後で気づきましたがauto_write=TRUEにしておけば大丈夫そうでした↩
原因を見つけてこの1行を追加するのに2時間くらいかかりました↩

サボタージュ禁止のおさぼり日記

NagoyaStat#8 アヒル本5章

NagoyaStat

ユーティリティの作成

Stanのコンパイル

予測値の分布

お礼