NIMBLEで並列化するための関数

MCMCをRで行うときに，たぶん一番良く使われるのはStanだと思います．
しかし，前の記事でも書きましたが，Stanはアルゴリズムの都合上，離散値パラメータの推定を行えません...
生態学の研究では（他の分野はどうなんだろう）個体数などの離散値パラメータが頻繁に登場するので，一々和を取って消すのはめんどくさいです．

なので，生態学では，シンプルに離散値パラメータを扱えるBUGS言語もよく使われます．
僕は学部時代からBUGS言語で書けるOpenBUGSを使い始めました．謎のエラー文と必死に戦ったのが懐かしいです．今はJAGSが主流だと思います．

これらに比べて，少しマイナーなNIMBLEというBUGS言語で書けるMCMCサンプラーがあります．
NIMBLEは，自作関数・確率分布が簡単に書けたり，推定に粒子フィルタが使えたりとかなり変態仕様です．しかも，WAICを一応自動で出してくれます．僕は研究で，ノンパラベイズ（ディリクレ過程関係）を使用したモデルを作っているので非常に重宝しています．

ただ，このNIMBLEは，StanやJAGSでは引数をちょっと変えるだけで使える並列化がデフォルトで使えません．
これは非常に不便で，自分で並列化するコードを書かなくてはなりません．

そこでこの記事では，NIMBLEを並列化するコードを紹介したいと思います．コードはこのサイトを参考に（ほぼ丸写し）しました．

並列化関数

いくつか方法があるようですが，foreachパッケージを利用します．

nimbleMCMC_para <- function(code, data, constants, inits, params, ni, nb, nt, nc){
  cl <- makeCluster(nc)
  registerDoParallel(cl)
  
  seeds <- 1:nc
  out0 <- foreach(x = seeds, .packages="nimble") %dopar% {
    set.seed(x)
    nimbleMCMC(code, data=data, constants=constants, inits=inits, monitors=params, niter = ni, nburnin = nb, thin = nt,nchains = 1, samplesAsCodaMCMC = TRUE) 
  }
  stopCluster(cl)
  out <- coda::mcmc.list(out0)
  return(out)
}

引数は基本的にnimbleMCMCと同じになっているはずです．チェイン毎に並列化するのでncの数だけコアが必要です．
nc = 8なのにコアが2個しかないみたいな状態でどうなるかは試してないのでわかりません...

戻り値は，mcmcサンプルが出てくるようにしてますが，この辺は好みによると思います．

NIMBLEはとにかく日本語の説明が少ないので，このブログで少しづつNIMBLEの魅力と使い方を紹介できたらと思います．