nimbleHMCを使ってみる

R
NIMBLE
作者

伊東宏樹

公開

2024年1月15日

NIMBLEでは、nimbleHMCパッケージにより、Stanで使われているNo-U-Turn Hamiltonian Monte Carloサンプラーを利用できます。ということで、試してみました。

モデル

例に、Eight schoolsのモデルを使います。

library(nimble)
## nimble version 1.0.1 is loaded.
## For more information on NIMBLE and a User Manual,
## please visit https://R-nimble.org.
## 
## Note for advanced users who have written their own MCMC samplers:
##   As of version 0.13.0, NIMBLE's protocol for handling posterior
##   predictive nodes has changed in a way that could affect user-defined
##   samplers in some situations. Please see Section 15.5.1 of the User Manual.
## 
##  次のパッケージを付け加えます: 'nimble'
##  以下のオブジェクトは 'package:stats' からマスクされています:
## 
##     simulate
library(nimbleHMC)
library(coda)

set.seed(123)

J <- 8
y <- c(28,  8, -3,  7, -1,  1, 18, 12)
sigma <- c(15, 10, 16, 11,  9, 11, 10, 18)

NIMBLEコード

NIMBLEのモデルコードです。

code <- nimbleCode({
  for (j in 1:J) {
    theta[j] <- mu + tau * eta[j]
    y[j] ~ dnorm(theta[j], sd = sigma[j])
    eta[j] ~ dnorm(0, sd = 1)
  }
  mu ~ dnorm(0, sd = 100)
  tau ~ dunif(0, 100)
})

初期値を生成しておきます。

inits <- list(eta = runif(J, -2, 2),
              mu = runif(1, -2, 2),
              tau = runif(1, 0, 2))

通常のMCMC

通常のMCMCによるパラメータ推定です。まず、nimbleModel関数でモデルを作成します。

model <- nimbleModel(code,
                     constants = list(J = J),
                     data = list(y = y, sigma = sigma),
                     inits = inits)
## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions

モデルをグラフで表示してみます。

model$plotGraph()

buildMCMC関数でMCMCオブジェクトを作成して、モデルとMCMCオブジェクトをコンパイルし、runMCMC関数でモデルをあてはめます。

mcmc <- buildMCMC(model,
                  monitors = c("theta", "mu", "tau"))
## ===== Monitors =====
## thin = 1: mu, tau, theta
## ===== Samplers =====
## RW sampler (1)
##   - tau
## conjugate sampler (9)
##   - eta[]  (8 elements)
##   - mu
compileNimble(model)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Derived CmodelBaseClass created by buildModelInterface for model code_MID_1
comp_mcmc <- compileNimble(mcmc, project = model)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
samp <- runMCMC(comp_mcmc,
                niter = 4000,
                nburnin = 2000,
                nchains = 3,
                samplesAsCodaMCMC = TRUE)
## running chain 1...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 2...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 3...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

収束診断

R-hatをチェックします。

gelman.diag(samp)
## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## mu             1.00       1.00
## tau            1.01       1.04
## theta[1]       1.00       1.00
## theta[2]       1.00       1.00
## theta[3]       1.00       1.01
## theta[4]       1.00       1.00
## theta[5]       1.01       1.02
## theta[6]       1.00       1.00
## theta[7]       1.00       1.00
## theta[8]       1.00       1.01
## 
## Multivariate psrf
## 
## 1.01

収束しているようです。

結果

推定結果です。

summary(samp)
## 
## Iterations = 1:2000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean    SD Naive SE Time-series SE
## mu        7.986 5.169  0.06673        0.10679
## tau       6.722 5.713  0.07376        0.26132
## theta[1] 11.591 8.515  0.10993        0.18681
## theta[2]  7.887 6.277  0.08104        0.08104
## theta[3]  5.976 7.936  0.10245        0.12532
## theta[4]  7.614 6.799  0.08777        0.08948
## theta[5]  5.128 6.458  0.08337        0.11287
## theta[6]  6.089 6.797  0.08775        0.09335
## theta[7] 10.766 6.924  0.08939        0.13917
## theta[8]  8.561 8.003  0.10332        0.11634
## 
## 2. Quantiles for each variable:
## 
##              2.5%   25%    50%    75% 97.5%
## mu        -2.0165 4.599  7.952 11.300 18.36
## tau        0.2282 2.569  5.339  9.363 21.25
## theta[1]  -2.2260 6.048 10.284 16.021 32.24
## theta[2]  -4.7943 3.948  7.934 11.729 20.63
## theta[3] -11.8518 1.848  6.511 10.785 20.71
## theta[4]  -6.1001 3.625  7.624 11.620 21.54
## theta[5]  -9.0236 1.269  5.720  9.514 16.62
## theta[6]  -8.5943 1.993  6.478 10.449 18.91
## theta[7]  -1.3919 6.168 10.062 14.697 26.30
## theta[8]  -7.2873 3.980  8.326 12.821 25.99

HMC

つづいて、HMCでパラメータ推定をしてみます。buildMCMCのかわりにbuildHMCを使います。引数にbuildDerivs = TRUEを指定します。

model_hmc <- nimbleModel(code,
                         constants = list(J = J),
                         data = list(y = y, sigma = sigma),
                         inits = inits,
                         buildDerivs = TRUE)
## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions
hmc <- buildHMC(model_hmc,
                monitors = c("theta", "mu", "tau"))
## ===== Monitors =====
## thin = 1: mu, tau, theta
## ===== Samplers =====
## NUTS sampler (1)
##   - eta[1], eta[2], eta[3], eta[4], eta[5], eta[6], eta[7], eta[8], mu, tau
compileNimble(model_hmc)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Derived CmodelBaseClass created by buildModelInterface for model code_MID_2
comp_hmc <- compileNimble(hmc, project = model_hmc)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
samp_hmc <- runMCMC(comp_hmc,
                    niter = 4000,
                    nburnin = 2000,
                    nchains = 3,
                    samplesAsCodaMCMC = TRUE)
## running chain 1...
##   [Note] NUTS sampler (nodes: eta[1], eta[2], eta[3], eta[4], eta[5], eta[6], eta[7], eta[8], mu, tau) is using 2000 warmup iterations.
##          Since `warmupMode` is 'default' and `nburnin` > 0,
##          the number of warmup iterations is equal to `nburnin`.
##          The burnin samples will be discarded, and all samples returned will be post-warmup.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 2...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 3...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

R-hat

R-hatをチェックします。

gelman.diag(samp_hmc)
## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## mu             1.00       1.00
## tau            1.01       1.01
## theta[1]       1.00       1.00
## theta[2]       1.00       1.00
## theta[3]       1.00       1.00
## theta[4]       1.00       1.00
## theta[5]       1.00       1.00
## theta[6]       1.00       1.00
## theta[7]       1.00       1.00
## theta[8]       1.00       1.00
## 
## Multivariate psrf
## 
## 1

こちらも収束しているようです。

結果

推定結果です。

summary(samp_hmc)
## 
## Iterations = 1:2000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean    SD Naive SE Time-series SE
## mu        7.954 5.013  0.06471        0.08868
## tau       6.537 5.770  0.07449        0.12433
## theta[1] 11.301 8.392  0.10835        0.12842
## theta[2]  7.906 6.263  0.08086        0.07054
## theta[3]  6.064 7.805  0.10076        0.10479
## theta[4]  7.710 6.619  0.08545        0.08241
## theta[5]  5.128 6.468  0.08350        0.07864
## theta[6]  5.989 6.693  0.08640        0.08644
## theta[7] 10.475 6.707  0.08658        0.09779
## theta[8]  8.331 7.871  0.10161        0.10034
## 
## 2. Quantiles for each variable:
## 
##              2.5%   25%    50%    75% 97.5%
## mu        -1.9213 4.668  7.924 11.079 18.17
## tau        0.2160 2.335  5.134  9.118 20.46
## theta[1]  -2.3830 5.775 10.153 15.522 31.96
## theta[2]  -4.7956 3.992  7.841 11.838 20.39
## theta[3] -11.4563 1.979  6.568 10.771 20.50
## theta[4]  -5.3690 3.697  7.676 11.681 21.35
## theta[5]  -9.1426 1.402  5.716  9.414 16.75
## theta[6]  -8.8994 2.138  6.329 10.297 18.39
## theta[7]  -0.8679 6.015  9.855 14.182 25.81
## theta[8]  -7.9434 3.861  8.059 12.447 25.34

nimbleHMC関数

nimbleHMC関数を使うと、nimbleMCMC関数と同様により簡単に実行できます。

samp_hmc2 <- nimbleHMC(code,
                       constants = list(J = J),
                       data = list(y = y, sigma = sigma),
                       inits = inits,
                       monitors = c("theta", "mu", "tau"),
                       niter = 4000,
                       nburnin = 2000,
                       nchains = 3,
                       samplesAsCodaMCMC = TRUE)
## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions
## Checking model calculations
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## running chain 1...
##   [Note] NUTS sampler (nodes: eta[1], eta[2], eta[3], eta[4], eta[5], eta[6], eta[7], eta[8], mu, tau) is using 2000 warmup iterations.
##          Since `warmupMode` is 'default' and `nburnin` > 0,
##          the number of warmup iterations is equal to `nburnin`.
##          The burnin samples will be discarded, and all samples returned will be post-warmup.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 2...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## running chain 3...
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

結果です。

summary(samp_hmc2)
## 
## Iterations = 1:2000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean    SD Naive SE Time-series SE
## mu        7.918 5.112  0.06599        0.13027
## tau       6.668 6.550  0.08456        0.35289
## theta[1] 10.960 8.414  0.10862        0.15904
## theta[2]  7.723 6.303  0.08137        0.10024
## theta[3]  6.192 7.674  0.09907        0.10527
## theta[4]  7.614 6.588  0.08505        0.07617
## theta[5]  5.118 6.214  0.08022        0.08486
## theta[6]  5.975 6.669  0.08609        0.08253
## theta[7] 10.638 6.837  0.08827        0.09298
## theta[8]  8.484 7.694  0.09934        0.14323
## 
## 2. Quantiles for each variable:
## 
##              2.5%   25%   50%    75% 97.5%
## mu        -2.0263 4.731 7.884 11.052 18.54
## tau        0.2248 2.343 5.198  9.022 21.05
## theta[1]  -3.3297 5.596 9.891 15.184 31.42
## theta[2]  -4.7184 3.879 7.757 11.543 20.59
## theta[3] -10.8754 2.105 6.659 10.829 20.59
## theta[4]  -5.7622 3.623 7.753 11.616 20.87
## theta[5]  -8.7050 1.409 5.420  9.389 16.40
## theta[6]  -9.0434 2.153 6.290 10.224 17.95
## theta[7]  -1.2866 6.063 9.944 14.563 25.92
## theta[8]  -6.6126 3.963 8.185 12.699 25.72