{cmdstanr}: Stan高速コンパイル、{rstan}代替としても有用

前回、最近mac版のRパッケージ、とくにコンパイル必要パッケージのインストールが困難になっているという記事を書いた。R3.6ではひどい状況だったが、R4.0.2の現在、CRANの手順に従えば、おおよそ問題はなくなっている。残る大きな問題は、rstanの動作不良だ。手元のR4.0.2では、rstanのインストール自体はできても、モデルのコンパイルができなくなっている。Macだけでなく、最近Windowsでもrstanの問題が出ている模様。

対策方法を調べていったところ、cmdstanrというコマンドラインベース（コンソール？macならTerminal）のRパッケージに辿り着いた。rstanのようにいちいちRcppでコンパイルしないので動作が速いし、またバージョンアップの際の依存環境との相性問題も生じにくいようだ。mac, linux, windowsのいずれでも使用可能なようだ。コマンドラインベースとは言うものの、Rから全ての作業を完結できることが分かったので、使い勝手も悪くない。

cmdstanrの公式ページが分かりやすいので、とくに説明も要らなさそうだが、実行サンプルを載せておく。

# インストール：　まず以下をR上で実行（公式ページの記述そのまま）

install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

# または、開発版のインストール（上記ではコンパイル失敗するモデルがあったが、そういう場合でも開発版では成功している2020.09.17追記 → 2021.09.30現在、通常インストールの最新版で十分にうまく行くようになった模様）

# おそらくもう不要→ # devtools::install_github("stan-dev/cmdstanr")

# 次に、以下でcmdstanrを呼び出し、

require(cmdstanr)

# R上で、下記のインストール関数を実行（coresの数値は、PCのコアをいくつ使用するかの指定、そこまで時間は掛からない）

install_cmdstan(cores = 2)

cmdstan_path() # インストールされるディレクトリを示す場合

# 2回目以降の使用で謎エラーが出てモデルがコンパイルできない場合、いったんcmdstanrを削除しインストールし直すことで使えている。いい対処法ではないだろうけれど、今のところ他にうまい手が見つからない。

####################################

# 計算実行：まず、モデル式を記述（公式サンプルではファイルからの読み込み）、

stan_program <- "

data {

  int<lower=0> N;

  int<lower=0,upper=1> y[N];

}

parameters {

  real<lower=0,upper=1> theta;

}

model {

  y ~ bernoulli(theta);

}

"

modfile <- write_stan_file(stan_program)

mod <- cmdstan_model(modfile) # 一時ファイルへ書き出し、この方が運用しやすいと思う

# あとはrstanと同様の計算実行（公式の計算例、そのまま）

data_list <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))

fit <- mod$sample(

  data = data_list,

  seed = 123, # 乱数の種

  chains = 4,

  parallel_chains = 4, # 並列化するchainの数

  refresh = 500 # コンソール上に経過を表示するiterationの更新間隔

)

# samplingの詳細設定の解説も充実

## 指定可能な代表的なもの

# iter_warmup: burn-in数

# iter_sampling: butn-in後のsampling数

# thin: sampleを保存する間隔

# init: 初期値、init=function() list(theta=runif(1))のように

# 計算結果はtibble形式になっている

fit$summary() # 計算結果の要約、tibble形式

fit$sampler_diagnostics() # 収束診断

fit$draws() # sample値（次元：iteration, chain, parameter）、

# "draws_array"形式だが、matrix(fit$draws()[,,1])あるいはmatrix(fit$draws("theta"))のようにparameter名で指定し、parameter毎に取り出し可能

# しかし、以下によってrstanの形式に変換することも可能

stanfit <- rstan::read_stan_csv(fit$output_files())

# rstanの形式にできた方が都合が良いだろうが、rstanのインストールが全くできない場合はtibbleの方がありがたい場合もあるだろう。

# 推定結果を保存する場合、以下のように.RDSで保存するのが推奨されている。

fit$save_object(file = "fit.RDS")

# 再読み込みの場合

fit2 <- readRDS("fit.RDS")

#（2021.12.31補足）個々のchainに１つのCPU threadを割り当てる並列化は上記のようにparallel_chainsの指定をするだけなので容易だが、chainあたり複数のthreadを割り当てるには、モデルコード自体にも工夫が要るようだ。GPUの利用なども可能なようだが、いずれにしてもモデル構造自体の制約が大きく、またN数が非常に大きいなどのビッグデータには効果が大きいが、生態学でよく用いるような複雑なプロセスモデルでは恩恵が少なそうである。https://mc-stan.org/docs/2_28/cmdstan-guide/parallelization.html