NIMBLEでMCMC

自己紹介

名前: 伊東宏樹

勤務先: 森林総合研究所北海道支所

共訳書:

NIMBLE

BUGS上位互換の言語で統計モデルを記述し、パラメータ推定ができるRパッケージ
モデルをC++に変換してコンパイルして実行する
2023年5月31日にVersion 1.0.0リリース（現在は1.0.1）
ウェブサイト: https://r-nimble.org/

インストール

C++コンパイラを使用するので、開発環境が必要¹

NIMBLE本体は普通にCRANからインストール

install.packages("nimble")

実行例

このデータに線形回帰モデルをあてはめ

モデルコード

nimbleCode関数の中でBUGS言語でモデルを記述

code <- nimbleCode({
  for (n in 1:N) {
    mu[n] <- alpha + beta * X[n]
    Y[n] ~ dnorm(mu[n], tau)
  }
  alpha ~ dnorm(0, 100)
  beta ~ dnorm(0, 100)
  tau <- 1 / (sigma * sigma)
  sigma ~ dunif(0, 100)
})

NIMBLEのBUGS拡張

code2 <- nimbleCode({
  for (n in 1:N) {
    Y[n] ~ dnorm(alpha + beta * X[n], sd = sigma)
  }
  alpha ~ dnorm(0, 100)
  beta ~ dnorm(0, 100)
  sigma ~ dunif(0, 100)
})

引数に式を使用可能
dnormのばらつきの指定に標準偏差を使用可能（分散も可）

MCMC

nimbleMCMC関数で、マルコフ連鎖モンテカルロ (MCMC) により、モデルのあてはめとパラメータ推定を実行

samp <- nimbleMCMC(code = code2,
                   constants = list(N = N),
                   data = list(X = X, Y = Y),
                   niter = 6000, nburnin = 1000,
                   nchains = 3,
                   samplesAsCodaMCMC = TRUE) |>
   posterior::as_draws()

最後にposteriorパッケージのdrawsクラスのオブジェクトに変換

結果の要約

summary(samp, default_summary_measures(), "rhat")

# A tibble: 3 × 8
  variable   mean median     sd    mad      q5   q95  rhat
  <chr>     <num>  <num>  <num>  <num>   <num> <num> <num>
1 alpha    0.0886 0.0887 0.0995 0.0995 -0.0749 0.252  1.00
2 beta     0.660  0.672  0.114  0.113   0.457  0.825  1.00
3 sigma    2.94   2.84   0.680  0.634   2.03   4.20   1.00

連鎖の軌跡のプロット

bayesplot::mcmc_trace(samp)

密度のプロット

bayesplot::mcmc_dens(samp)

その他の機能

WAIC (Widely Applicable Information Criterion)
空間統計モデル（CAR: Conditional AutoRegression）
逐次モンテカルロ（粒子フィルタ, nimbleSMCパッケージ）

など、機能豊富（詳細は公式マニュアルを）

🥳