Working through quick-start example in nimble manual

The manual gives essentially no introduction to what appears to be a classic BUGS example model for stochastically failing pumps.

```
library(nimble)
pumpCode <- modelCode({
for (i in 1:N){
theta[i] ~ dgamma(alpha,beta)
lambda[i] <- theta[i]*t[i]
x[i] ~ dpois(lambda[i])
}
alpha ~ dexp(1.0)
beta ~ dgamma(0.1,1.0)
})
```

```
pumpConsts <- list(N = 10,
t = c(94.3, 15.7, 62.9, 126, 5.24,
31.4, 1.05, 1.05, 2.1, 10.5))
pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
```

```
pumpInits <- list(alpha = 1,
beta = 1,
theta = rep(0.1, pumpConsts$N))
```

```
pump <- nimbleModel(code = pumpCode,
name = 'pump',
constants = pumpConsts,
data = pumpData,
inits = pumpInits)
pump$getNodeNames()
```

```
[1] "alpha" "beta" "lifted_d1_over_beta"
[4] "theta[1]" "theta[2]" "theta[3]"
[7] "theta[4]" "theta[5]" "theta[6]"
[10] "theta[7]" "theta[8]" "theta[9]"
[13] "theta[10]" "lambda[1]" "lambda[2]"
[16] "lambda[3]" "lambda[4]" "lambda[5]"
[19] "lambda[6]" "lambda[7]" "lambda[8]"
[22] "lambda[9]" "lambda[10]" "x[1]"
[25] "x[2]" "x[3]" "x[4]"
[28] "x[5]" "x[6]" "x[7]"
[31] "x[8]" "x[9]" "x[10]"
```

Note that we can see `theta`

has our initial conditions, while `lambda`

has not yet been initialized:

`pump$theta`

` [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1`

`pump$lambda`

` [1] NA NA NA NA NA NA NA NA NA NA`

Hmm, initially we cannot simulate `theta`

values though (or rather, we just get NaNs and warnings if we do). At the moment I’m not clear on why, though seems to be due to the lifted node:

```
simulate(pump, 'theta')
pump$theta
```

` [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN`

`pump$lifted_d1_over_beta`

`[1] NA`

If we calculate the log probability density of the determinstic dependencies of alpha and beta nodes (i.e. the lifted node) then we’re okay:

```
set.seed(0) ## This makes the simulations here reproducible
calculate(pump, pump$getDependencies(c('alpha', 'beta'), determOnly = TRUE))
```

`[1] 0`

`simulate(pump, 'theta')`

`pump$theta`

```
[1] 1.79180692 0.29592523 0.08369014 0.83617765 1.22254365 1.15835525
[7] 0.99001994 0.30737332 0.09461909 0.15720154
```

We still need to initialize lambda, e.g. by calculating the probability density on those nodes:

`calculate(pump, 'lambda')`

`[1] 0`

`pump$lambda`

```
[1] 168.9673926 4.6460261 5.2641096 105.3583839 6.4061287
[6] 36.3723548 1.0395209 0.3227420 0.1987001 1.6506161
```

though not entirely clear to me why the guide prefers to do this as the dependencies of theta (which clearly include lambda, but also other things). Also not clear if these `calculate`

steps are necessary to proceed with the `MCMCspec`

and `buildMCMC`

, or compile steps. Let’s reset the model^{1} and find out:

```
pump <- nimbleModel(code = pumpCode,
name = 'pump',
constants = pumpConsts,
data = pumpData,
inits = pumpInits)
pump$theta
```

` [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1`

`pump$lambda`

` [1] NA NA NA NA NA NA NA NA NA NA`

Good, we’re reset. Now we try:

```
Cpump <- compileNimble(pump)
pumpSpec <- MCMCspec(pump)
pumpSpec$addMonitors(c('alpha', 'beta', 'theta'))
```

`thin = 1: alpha, beta, theta`

```
pumpMCMC <- buildMCMC(pumpSpec)
CpumpMCMC <- compileNimble(pumpMCMC, project = pump)
```

`CpumpMCMC(1000)`

`NULL`

`samples <- as.matrix(nfVar(CpumpMCMC, 'mvSamples'))`

```
plot(samples[ , 'alpha'], type = 'l', xlab = 'iteration',
ylab = expression(alpha))
plot(samples[ , 'beta'], type = 'l', xlab = 'iteration',
ylab = expression(beta))
plot(samples[ , 'alpha'], samples[ , 'beta'], xlab = expression(alpha),
ylab = expression(beta))
```

Note the poor mixing (which is improved by the block sampler, as shown in the manual).

Not completely certain that this destroys anything connected to the object as C pointers from before, but seems like it should.↩