NIMBLEでノンパラベイズを試したメモ

NIMBLEというRのライブラリがあります。BUGS言語風の文法でC++にコンパイルされるタイプの確率的プログラミング言語です。実装されている推定のアルゴリズムはここに書いてあります。MCMCの他にも以下のようなアルゴリズムがデフォルトで実装されており、実行速度もかなり速いです。

• particle filter（Sequential Monte Carlo, SMC）
• 空間モデルのCARモデル（Gaussian Markov Random Fieldsと等価）
• ノンパラベイズのChinese Restaurant Process（CRP）とStick-breaking Process（SBP）

ここでは、ノンパラベイズのCRPとSBPを試してみます。参考にしたのは公式ドキュメントです。

データは以下の以前の記事と同じで、モデルも似たようなものです。

statmodeling.hatenablog.com

はじめにCRP版、次にSBP版を紹介します。

Chinese Restaurant Process（CRP）

library(nimble)

code <- nimbleCode({
  group[1:N] ~ dCRP(alpha, size=N)
  alpha ~ dgamma(1, 1)
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

set.seed(1)
load('data-nonpara-bayes.RData')
constants <- list(N = data$n, C = data$C)
data <- list(Y = data$velocity)
inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              b = 1, alpha = 1)

fit <- nimbleMCMC(code = code, constants = constants,
                  data = data, inits = inits,
                  nchains = 4, niter = 10000, nburnin = 2000,
                  summary = TRUE, samplesAsCodaMCMC = TRUE,
                  monitors = c('group', 'mu_mix', 's2_mix' , 'alpha'))

library(coda)
pdf(file='fit-traceplot-CRP.pdf')
traceplot(fit$samples)
dev.off()

• 4行目の group[1:N] ~ dCRP(alpha, size=N) でCRPに従うサンプリングができます。専用のアルゴリズムになっているとのことで、BUGSやStanに比べると動作がだいぶ高速です。
• 6行目の C はクラスター数で、Nよりはだいぶ小さく、ただし、クラスターとクラスターがつぶれて合体しないように少し大きめにとる必要があります（nimbleを実行する際にその旨のwarningが表示されます）。データ（data-nonpara-bayes.RData）に含まれていますが、ここでは20に設定しています。
• 20行目　Stanと異なり、全てのパラメータの初期値をもれなくきちんと与えておかないと、尤度が無限やNAになってRごと落ちたりしますので注意。
• 28行目 samplesAsCodaMCMC = TRUEを設定しておくと{coda}ライブラリで扱いやすい形のサンプリング結果となり、31～34行目のようにしてtraceplotを出力することができます。

その他にも、コンパイルとサンプリングの実行を別にできたりしますが、ドキュメントが若干分かりにくいです。

Stick-breaking Process（SBP）

library(nimble)

code <- nimbleCode({
  for(c in 1:(C-1)) {
    v[c] ~ dbeta(1, alpha)
  }
  alpha ~ dgamma(1, 1)
  w[1:C] <- stick_breaking(v[1:(C-1)])
  for(c in 1:C) {
    mu_mix[c] ~ dnorm(20, 100)
    s2_mix[c] ~ dinvgamma(5, b)
  }
  b ~ dgamma(0.01, 0.01)
  for(n in 1:N) {
    group[n] ~ dcat(w[1:C])
    Y[n] ~ dnorm(mu_mix[group[n]], var=s2_mix[group[n]])
  }
})

(...CRP版と同じ...)

inits <- list(mu_mix = seq(0, 20, len = constants$C),
              s2_mix = rep(4, constants$C),
              group = sample(1:10, size = constants$N, replace = TRUE),
              v = rbeta(constants$C, 1, 1),
              b = 1, alpha = 1)

(...CRP版と同じ...)

CRP版とほとんど同じです。

どちらの場合においても、サンプリング結果からクラスター数や密度を算出する関数を自作する必要があります。またモデルのコードにバグが含まれていると、コンパイルエラーではなくRごと落ちたりするのは若いライブラリのご愛敬です。多少不便でもMCMCの性能はかなり良いと思うので積極的に使いたいと思います。

Enjoy!