A quick first exploration of NIMBLE and some questions.

```
library("nimble")
library("sde")
```

Let’s simulate from a simple OU process: \(dX = \alpha (\theta - X) dt + \sigma dB_t\)

```
set.seed(123)
d <- expression(0.5 * (10-x))
s <- expression(1)
data <- as.data.frame(sde.sim(X0=6,drift=d, sigma=s, T=100, N=400))
```

`## sigma.x not provided, attempting symbolic derivation.`

i.e. \(\alpha = 0.5\), \(\theta = 10\), \(\sigma=1\), starting at \(X_0 = 6\) and running for 100 time units with a dense sampling of 400 points.

Le’t now estimate a Ricker model based upon (set aside closed-form solutions to this estimate for the moment, since we’re investigating MCMC behavior here).

```
code <- modelCode({
K ~ dunif(0.01, 40.0)
r ~ dunif(0.01, 20.0)
sigma ~ dunif(1e-6, 100)
iQ <- 1 / (sigma * sigma)
x[1] ~ dunif(0, 10)
for(t in 1:(N-1)){
mu[t] <- log(x[t]) + r * (1 - x[t]/K)
x[t+1] ~ dlnorm(mu[t], iQ)
}
})
constants <- list(N = length(data$x))
inits <- list(K = 6, r = 1, sigma = 1)
Rmodel <- nimbleModel(code=code, constants=constants, data=data, inits=inits)
```

NIMBLE certainly makes for a nice syntax so far. Here we go now: create MCMC specification and algorithm

```
mcmcspec <- MCMCspec(Rmodel)
Rmcmc <- buildMCMC(mcmcspec)
```

Note that we can also query some details regarding our specification (set by default)

`mcmcspec$getSamplers()`

```
## [1] RW sampler; targetNode: K, adaptive: TRUE, adaptInterval: 200, scale: 1
## [2] RW sampler; targetNode: r, adaptive: TRUE, adaptInterval: 200, scale: 1
## [3] RW sampler; targetNode: sigma, adaptive: TRUE, adaptInterval: 200, scale: 1
```

`mcmcspec$getMonitors()`

`## thin = 1: K, r, sigma, x`

Now we’re ready to compile model and MCMC algorithm

```
Cmodel <- compileNimble(Rmodel)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
```

Note we could have specified the `Rmodel`

as the “project” (as shown in the example from the Nimble website), but this is more explicit. Rather convenient way to add to an existing model in this manner.

And Now we can execute the MCMC algorithm in blazing fast C++ and then extract the samples

`Cmcmc(10000)`

`## NULL`

`samples <- as.data.frame(as.matrix(nfVar(Cmcmc, 'mvSamples')))`

How do these estimates compare to theory:

`mean(samples$K)`

`## [1] 10.05681`

`mean(samples$r)`

`## [1] 0.180207`

Some quick impressions:

Strange that

`Rmodel`

call has to be repeated before we can set up a custom MCMC (nimble docs). How/when was this object altered since it was defined in the above code? Seems like this could be problematic for interpreting / reproducing results?What’s going on with

`getSamplers()`

and`getMonitors()`

? Guessing these are in there just to show us what the defaults will be for our model?why do we assign

`Cmodel`

if it seems we don’t use it? (the compilation needs to be done but isn’t explicitly passed to the next step). Seems we can use`Cmodel`

instead of`Rmodel`

in the`Cmcmc <- compileNimble(Rmcmc, project = Cmodel)`

, which makes the dependency more explicit, at least that notation is more explicit. Seems like it should be possiple to omit the first`compileNimble()`

and have the second call the`compileNimble`

automatically if it gets an object whose class is that of`Rmodel`

instead?Repeated calls to

`Cmcmc`

seem not to give the same results. Are we adding additional mcmc steps by doing this?Thinking an

`as.data.frame`

would be nicer than`as.matrix`

in the`nfVar`

`mvSamples`

coercion.Don’t understand what

`simulate`

does (or why it always returns NULL?).