R: NIMBLE の使い方 (その1)

こちらのページと同じことを行いました。
Introduction to NIMBLE
冬前に、57匹の動物を放し、冬が開けたら、19匹の動物を捕まえたという状況の分析です。

ライブラリーのインストール

install.packages('nimble')
install.packages('MCMCvis')

プログラム

animals.R

#
library(nimble)
library(MCMCvis)
#
computeLifespan <- nimbleFunction(
    run = function(theta = double(0)) { # type declarations
        ans <- -1/log(theta)
        return(ans)
        returnType(double(0))  # return type declaration
    } )
#
computeLifespan(0.8)
# nCcomputeLifespan <- compileNimble(computeLifespan)
# CcomputeLifespan(0.8)
#
model <- nimbleCode({
  # likelihood
  survived ~ dbinom(theta, released)
  # prior
  theta ~ dunif(0, 1)
  # derived quantity
  lifespan <- computeLifespan(theta)
})
#
my.data <- list(survived = 19, released = 57)
parameters.to.save <- c("theta", "lifespan")
initial.values <- function() list(theta = runif(1,0,1))
n.iter <- 5000
n.burnin <- 1000
n.chains <- 2
mcmc.output <- nimbleMCMC(code = model,
                          data = my.data,
                          inits = initial.values,
                          monitors = parameters.to.save,
                          niter = n.iter,
                          nburnin = n.burnin,
                          nchains = n.chains)
MCMCsummary(object = mcmc.output, round = 2)
#

実行結果

$ Rscript animals.R 
nimble version 0.12.2 is loaded.
For more information on NIMBLE and a User Manual,
please visit https://R-nimble.org.

 次のパッケージを付け加えます: ‘nimble’ 

 以下のオブジェクトは ‘package:stats’ からマスクされています:

    simulate

[1] 4.48142
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 ...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
Running chain 2 ...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
         mean   sd 2.5%  50% 97.5% Rhat n.eff
lifespan 0.93 0.16 0.67 0.92  1.28    1  1799
theta    0.34 0.06 0.22 0.34  0.46    1  1818

次のバージョンで確認しました。

$ Rscript --version
Rscript (R) version 4.2.1 (2022-06-23)