Monday_NIMBLE.Rmd

---
title: "Monday NIMBLE"
author: "Mike Li, Morgan Kain"
date:  "`r format(Sys.time(), '%H:%M %d %B %Y')`"
output: html_document
---

```{r setup}
library(knitr)
options(mc.cores = parallel::detectCores())
#opts_chunk$set(cache = TRUE)
```

## NIMBLE: Numerical Inference for statistical Models for Bayesian and Likelihood Estimation

NIMBLE is built in R but compiles your models and algorithms using C++ for speed <br />
NIMBLE is most commonly used for MCMC but can also be used to implement a series of other algorithms (e.g. particle filtering, MCEM) <br />
<br />
1. A system for writing statistical models flexibly, which is an extension of the BUGS language <br />
2. A library of algorithms such as MCMC. <br />
3. A language, called NIMBLE, embedded within and similar in style to R, for writing algorithms that operate on BUGS models. <br />
  
One of the most important concepts behind NIMBLE is to allow a combination of highlevel processing in R and low-level processing in compiled C++. <br />
<br />

##### Why NIMBLE?

1. Options (More customizable MCMC, ability to run JAGS models and STAN models, EM, particle filter) that leads to a more adaptable workflow <br />
2. User-defined functions and distributions – written as nimbleFunctions – can be used in model code.  <br />
3. Multiple parameterizations for distributions, similar to those in R, can be used. <br />
<br />
  e.g. normal distribution with BUGS parameter order: <br />
        x ~ dnorm(a + b * c, tau) <br />
       normal distribution with a named parameter: <br />
        y ~ dnorm(a + b * c, sd = sigma) <br />
<br />
4. Named parameters for distributions and functions, similar to R function calls, can be used. <br />
5. More flexible indexing of vector nodes within larger variables is allowed. For example one can place a multivariate normal vector arbitrarily within a higher-dimensional object, not just in the last index. <br />
6. More general constraints can be declared using dconstraint, which extends the concept of JAGS’ dinterval. <br />
<br />

#### Downloading, installing and loading NIMBLE

On Windows, you should download and install Rtools.exe available from http://cran. r-project.org/bin/windows/Rtools/.  <br />
On OS X, you should install Xcode.  <br />

After these are installed you can install NIMBLE in R using <br />
install.packages("nimble", repos = "http://r-nimble.org", type = "source") <br />

Please post about installation problems to the nimble-users Google group or email nimble.stats@gmail.com.

You will also need to download STAN using the following commands <br />
Sys.setenv(MAKEFLAGS = "-j4") <br />
install.packages("rstan", dependencies = TRUE) <br />
<br />

In total you will need the following pakages:

```{r loadpackage, echo=TRUE, message=FALSE, verbose = FALSE}
library("nimble")
library("R2jags")
library("ggplot2")
library("nimble")
library("rstan")
library("igraph")
library("parallel")
library("mcmcplots")
library("lattice")
library("coda")
library("reshape2")
```

<br />
<br />

#### Things to know about working with NIMBLE

<b> Programming in NIMBLE involves a fundamental distinction between: </b> <br />
  1. the steps for an algorithm that need to happen only once, at the beginning, such as inspecting the model <br />
  2. the steps that need to happen each time a function is called, such as MCMC iterations. <br />
  <br />
    When one writes a nimbleFunction, each of these parts can be provided separately. 
<br />

Multiple parameterizations for distributions, similar to those in R, can be used.
NIMBLE calls non-stochastic nodes “deterministic”, whereas BUGS calls them “logical”. 
NIMBLE uses “logical” in the way R does, to refer to boolean (TRUE/FALSE) variables. <br />
Alternative models can be defined from the same model code by using if-then-else statements that are evaluated when the model is defined.

1. NIMBLE extracts all the declarations in the BUGS code to create a model definition. <br />
2. From the model definition, NIMBLE builds a working model in R. This can be used to manipulate variables and operate the model from R. Operating the model includes calculating, simulating, or querying the log probability value of model nodes. <br />
3. From the working model, NIMBLE generates customized C++ code representing the model, compiles the C++, loads it back into R, and provides an R object that interfaces to it. We often call the uncompiled model the “R-model” and the compiled model the “C-model.” <br />  
<br />

### Presentation Outline
The general outline for this presentation follows along with the NIMBLE users manual <br />
http://r-nimble.org/documentation-2 <br />
However, the model(s) used here are written by us <br />

##### Part 1
[1.1](#1.1) Build a chain binomial model in JAGS. Conduct parameter estimation <br />
[1.2](#1.2) Translate the model into NIBLE. Conduct parameter estimation <br />
\ \ \ \ \ [1.2.1](#1.2.1) Model exploration/conversion <br />
\ \ \ \ \ [1.2.2](#1.2.2) Create a basic MCMC specification for the chain binomial, compile and run the MCMC <br />
\ \ \ \ \ [1.2.3](#1.2.3) Small MCMC specification adjustments (more on this in Part 3) <br />
[1.3](#1.3) Compare the JAGS and NIMBLE results (parameter estimates, uncertainty, convergence, efficiency) <br />

##### Part 2
[2.1](#2.1) Translate the model using a "hybrid approach" (STAN does not allow for discrete latent variables) <br />
\ \ \ \ \ [2.1.1](#1.4.1) Conduct parameter estimation using JAGS and NIMBLE <br />
\ \ \ \ \ [2.1.2](#1.4.2) Run the hybrid model in STAN and compare the results from JAGS, NIMBLE and STAN <br />
[2.2](#2.2) Compare the NIMBLE Chain Binomial and STAN hybrid model <br />

##### Part 3
[3.1](#3.1) Expolore more fine-tuned adjustments that can be made in NIMBLE <br />
\ \ \ \ \ [3.1.1](#3.1.1)  NIMBLE functions (e.g. allows for the implementation of custom samplers) <br />
      
##### Part 4     
[4.1](#4.1) NIMBLE extras: <br />
\ \ \ \ \ [4.1.1](#4.1.1) Create, compile and run a Monte Carlo Expectation Maximization (MCEM) algorithm, which illustrates some of the flexibility NIMBLE provides to combine R and NIMBLE. <br />
\ \ \ \ \ [4.1.2](#4.1.2) Implement particle filtering for the chain binomial <br />
      
##### Part 5

[5.1](#5.1) Misc NIMBLE notes (truncated distributions, lifted nodes, logProb, multiple instances of the same model)

<br />

### Part 1

##### <a name="1.1"> 1.1 Build a chain binomial model in JAGS </a>
First step is to construct the simulator from which we will obtain our data <br />

Note: It will be important to set your current working directory to "../stat744/notes/NIMBLE" <br />

Set parameters and load the Chain Binomial simulator <br />

```{r, cache.rebuild = TRUE}
beta <- 0.02
pop <- 100
effpropS <- 0.8
effpropI <- 0.2
reporting <- 0.5

s0 <- effpropS*pop
r0 <- 0
zerohack <- 0.001
numobs <- 12
nimtimevec <- c()
source("CBsimulatorMK.R")
```


```{r, cache.rebuild = TRUE}
sim <- simCB(beta = beta, pop = pop, effpropS = effpropS, effpropI = effpropI, 
             t0 = 1, numobs = numobs, reporting = reporting, seed = 5)
sim
```

Take a peek at what this model produces

```{r, echo=FALSE}
  ggplot(data.frame(sim), aes(time, S)) + geom_point(colour = "blue") + 
    geom_line(colour = "blue") +
    geom_line(data = data.frame(sim), aes(time, I), colour = "red") +
    geom_line(data = data.frame(sim), aes(time, R), colour = "green") +
    theme_bw()
```

Set up the required arguments to run the JAGS model

```{r}
data <- list(obs = sim$Iobs,
             pop = pop,
             numobs = nrow(sim),
             r0 = r0)

inits <- list(list(
  I = sim$I*1,
  effpropS = effpropS - 0.1,
  effpropI = effpropI - 0.15,
  beta = beta + .05,
  reporting = reporting),
list(
  I = sim$I*1 + 1,
  effpropS = effpropS - 0.1,
  effpropI = effpropI + 0.2,
  beta = beta + .1,
  reporting = reporting),
list(
  I = sim$I*1 + 2,
  effpropS = effpropS,
  effpropI = effpropI,
  beta = beta - .008,
  reporting = reporting)
)

params = c("beta",
           "effpropS",
           "effpropI",
           "reporting")

#rjags::set.factory("bugs::Conjugate", FALSE, type="sampler")
```

Create the model and examine the MCMC algorithms that JAGS will use to sample <br />

```{r}
cbjagsmodel <- jags.model(data = data,
               inits = inits,
               file = "CB.bug",
               n.chains = length(inits))

list.samplers(cbjagsmodel)
```

Run some chains (could use coda::coda.samples from cbjagsmodel but here we will just run jags()) <br />

```{r, eval = FALSE}
coda.samples(cbjagsmodel, c("beta", "effpropS", "effpropI", "reporting"), 1000)
```

Run through jags <br />

```{r}
jagstime <- system.time(cbjags <- jags(data = data,
               inits = inits,
               param = params,
               model.file = "CB.bug",
               n.iter = 11000,
               n.burnin = 1000,
               n.thin = 20,
               n.chains = length(inits)))
```

```{r}
cbjags
xyplot(as.mcmc(cbjags))
```

<br />
<br />

##### <a name="1.2"> 1.2 Build the NIMLE model </a>

```{r}
source('nimCBMK.R')
```

Set up the model. Here we need: Constants, Data, Initial Values, NIMBLE model object <br />

```{r}
nimCBdata <- list(obs = sim$Iobs)

nimCBcon <- list(numobs = numobs, pop = pop, r0 = r0)

nimCBinits <- list(I = sim$I,
                   effpropS = effpropS,
                   effpropI = effpropI,
                   beta = beta,
                   reporting = reporting,
                   s0 = s0)

nimtimevec[1] <- system.time(CBout <- nimbleModel(code = nimcode, 
                         name = 'CBout', 
                         constants = nimCBcon,
                         data = nimCBdata, 
                         inits = nimCBinits))[3]
```

<br />
<br />

##### <a name="1.2.1"> 1.2.1 Model exploration/conversion </a>

```{r}
CBout$getNodeNames()
```

You can examine a single specific node using: <br />

```{r}
CBout[['I[2]']]
```

```{r}
CBout$obs
```

nimbleModel does its best to initialize a model, but let’s say you want to re-initialize I. <br />

```{r}
simulate(CBout, 'I') # using the current beta -- if we update beta to a new value this will change
CBout$I
```

And take a look at the log-prob <br />

```{r}
calculate(CBout, "I")
CBout$logProb_I

I2lp <- CBout$nodes[['I[2]']]$calculate()
I2lp
```

or Calculate new log probabilities after updating I <br />

```{r}
CBout$obs

getLogProb(CBout, 'obs')

calculate(CBout, CBout$getDependencies(c("obs")))

CBout$logProb_obs
```

We can also look at the dependencies to make sense of what is going on <br />

```{r}
par(mfrow = c(1,1))
plot(CBout$graph)
```

```{r}
CBout$getDependencies(c("beta"))
```

```{r}
CBout$getDependencies(c("beta"), determOnly = TRUE)
```

We will return to nimbleFunctions later -- but here is an initial example nimbleFunction that will simulate multiple values for a designated set of nodes and calculate every part of the model that depends on them <br />

```{r}
simNodesMany <- nimbleFunction( setup = function(model, nodes) {
        mv <- modelValues(model)
        deps <- model$getDependencies(nodes)
        allNodes <- model$getNodeNames()
},
run = function(n = integer()) {
  resize(mv, n) 
  for(i in 1:n) {
            simulate(model, nodes)
            calculate(model, deps)
            copy(from = model, nodes = allNodes,
to = mv, rowTo = i, logProb = TRUE) }
})
```

```{r}
outI <- simNodesMany(CBout, "I")

outI$run(1)

outI$mv[["logProb_obs"]]
```

<br />
<br />

##### <a name="1.2.2"> 1.2.2 Compile the C++ model, build the MCMC and Run </a>

```{r}
nimtimevec[2] <- system.time(CBoutC <- compileNimble(CBout))[3]
```

Configure the MCMC with the default options (we will return to customizing this setup later) <br />

```{r, message=FALSE, verbose = FALSE, warning = FALSE, errors = FALSE}
nimtimevec[3] <- system.time(CBoutSpec <- configureMCMC(CBout, print = TRUE))[3]
```

Add chain monitors for the parameters of interest and add thinning <br />

```{r}
CBoutSpec$addMonitors(c("beta", "effpropS", "effpropI", "reporting"))

CBoutSpec$setThin(20)
```

Build the MCMC <br />

```{r}
nimtimevec[4] <- system.time(CBoutMCMC <- buildMCMC(CBoutSpec))[3]

nimtimevec[5] <- system.time(CBoutMCMC <- compileNimble(CBoutMCMC, project = CBout, resetFunctions = TRUE))[3]
```

```{r}
niter <- 11000

set.seed(0)
nimtimevec[6] <- system.time(CBoutMCMC$run(niter))[3]
  # mcmc$run(niter, reset = FALSE) can be used to add more iterations
```

Quick peek at time required. Below we will take a look at efficiency using a convenient coda function <br />

Gross time required <br />

```{r}
jagstime[3]
sum(nimtimevec[1:6], na.rm = TRUE)
nimtimevec[6]
```

Efficiency (Net time in a sense) <br />

```{r}
samples <- as.matrix(CBoutMCMC$mvSamples)

head(samples)

jags_eff <- effectiveSize(as.mcmc.list(as.mcmc(cbjags))) / nimtimevec[1]
nim_eff <- effectiveSize(as.mcmc.list(as.mcmc(samples))) / nimtimevec[6]

jags_eff
nim_eff
```


```{r, echo = FALSE}
  samples <- as.matrix(CBoutMCMC$mvSamples)
  par(mfrow = c(2, 2), mai = c(.5, .5, .3, .2))
  plot(samples[ , 'beta'], type = 'l', xlab = 'iteration', main = "beta")
  plot(samples[ , 'effpropS'], type = 'l', xlab = 'iteration', main = "effpropS")
  plot(samples[ , 'effpropI'], type = 'l', xlab = 'iteration', main = "effpropI")
  plot(samples[ , 'reporting'], type = 'l', xlab = 'iteration', main = "reporting")
```

Look at the correlation in the chains

```{r}
  acf(samples[, "beta"])
  acf(samples[, "reporting"])
  acf(samples[, "effpropS"])
  acf(samples[, "effpropI"])
```

<br />
<br />

##### <a name="1.2.3"> 1.2.3 Small MCMC specification adjustment </a>

A few undesirable results here... we can add a block sampler to decrease correlation <br />

Take a look at the samplers being used <br />

```{r}
CBoutSpec$getSamplers()
```
  
```{r}
CBoutSpec$addSampler(target = c('effpropS', 'effpropI'), type = 'RW_block',
                      control = list(adaptInterval = 10000))
```

```{r}
CBoutSpec$setThin(20)
```

```{r}
CBoutMCMC <- buildMCMC(CBoutSpec)
```

```{r}
CBoutMCMC <- compileNimble(CBoutMCMC, project  = CBout, resetFunctions = TRUE)
```

```{r}
nimtimevec[7] <- system.time(CBoutMCMC$run(11000))[3]
samplesNew <- as.matrix(CBoutMCMC$mvSamples)
```

Check for an imporvement

```{r}
  acf(samplesNew[, "effpropS"])
  acf(samplesNew[, "effpropI"])
  plot(samplesNew[ , 'effpropS'], type = 'l', xlab = 'iteration')
  plot(samplesNew[ , 'effpropI'], type = 'l', xlab = 'iteration')
```

```{r}
nim_eff2 <- effectiveSize(as.mcmc.list(as.mcmc(samplesNew))) / nimtimevec[7]

nim_eff
nim_eff2
```

NIMBLE allows for specification of samplers by parameter or node by node (NIMBLE included or user created) <br />

e.g. <br />

```{r, eval = FALSE}
CBoutSpec$removeSamplers(c("beta", "effpropS", "effpropI", "reporting"), print = TRUE)
CBoutSpec$addSampler("beta", type = "slice", print = TRUE)
CBoutSpec$addSampler("effpropS", type = "slice", print = TRUE)
CBoutSpec$addSampler("effpropI", type = "slice", print = TRUE)
CBoutSpec$addSampler("reporting", type = "slice", print = TRUE)
```

<br />
<br />

##### <a name="1.3"> 1.3 Compare the JAGS and NIMBLE results </a>

We can also compare the NIMBLE model simultaneously with the JAGS model using MCMCsuite() <br />

Be warned: running this code with "makePlot = TRUE" will produce about 6-8 graphs which will all pop up in separate windows! <br />

```{r, "Homework"}
source("simulateCB.R")
options(mc.cores = parallel::detectCores())
source('nimCB.R')

nimCBdata <- list(obs=sim$Iobs)
nimCBcon <- list(numobs=numobs,N=N,i0=i0)

nimCBinits <- list(I=sim$I,
                   effprop=effprop,
                   beta=beta,
                   reporting=reporting,
                   N0=N0
)
NimbleCB <- MCMCsuite(code=nimcode,
                   data=nimCBdata,
                   inits=nimCBinits,
                   constants=nimCBcon,
                   MCMCs=c("jags","nimble"),
                   monitors=c("beta","reporting","effprop"),
                   calculateEfficiency=TRUE,
                   niter=iterations,
                   makePlot=FALSE,
                   savePlot=FALSE,
                   setSeed = 5)

print(NimbleCB$timing)
print(NimbleCB$summary)

```

It is kind of scary considering Nimble is suppose to be this awesome super fast magic box. Before we learn Nimble, let's take a few minutes learning JAGS and samplers.


```{r, Jags Example}

rjags::set.factory("bugs::Conjugate", FALSE, type="sampler")

Jagsmod <- jags.model(file="CB2.bug",data=data,inits=inits)

list.samplers(Jagsmod)

slicetime <- system.time(JagsCB <- jags(data=data,
               inits=inits,
               param = params,
               model.file = "CB2.bug",
               n.iter = iterations,
               n.chains = length(inits),
               n.thin = 1,
               n.burnin = 2000
               ))

print(JagsCB)
slice_eff <- effectiveSize(as.mcmc(JagsCB))/slicetime[3]


rjags::set.factory("bugs::Conjugate", TRUE, type="sampler")

Jagsmod2 <- jags.model(file="CB2.bug",data=data,inits=inits)

list.samplers(Jagsmod2)

conjutime <- system.time(JagsCB2 <- jags(data=data,
               inits=inits,
               param = params,
               model.file = "CB2.bug",
               n.iter = iterations,
               n.chains = length(inits),
               n.thin = 1,
               n.burnin = 2000
               ))
print(JagsCB)
print(JagsCB2)
conju_eff <- effectiveSize(as.mcmc(JagsCB2))/conjutime[3]

slicetime
conjutime

slice_eff
conju_eff
```

Wow, that is pretty cool. Let's try to do it in Nimble. 

```{r, NimbleCB}
NimbleCB <- MCMCsuite(code=nimcode,
                   data=nimCBdata,
                   inits=nimCBinits,
                   constants=nimCBcon,
                   MCMCs=c("jags","nimble","nimble_slice"),
                   monitors=c("beta","reporting","effprop"),
                   calculateEfficiency=TRUE,
                   niter=iterations,
                   makePlot=FALSE,
                   savePlot=FALSE,
                   setSeed = 5)

print(NimbleCB$timing)
print(NimbleCB$summary)

```

Why is default Nimble underperforming? (vs Jags) Take my word, nimble-jags is NOT using conjugate samplers. Nimble slice turns all nodes to slice samplers, thus, nimble-jags and nimble slice is using the same sampler. As you can see Nimble slice is more efficient than JAGS.

```{r, "Investigating default Nimble"}
mod <- nimbleModel(code = nimcode, data=nimCBdata, inits=nimCBinits, constants=nimCBcon, name = "mod")
temp <- compileNimble(mod) ## need to compile it once before recompiling changes
Cmod <- configureMCMC(mod,print=TRUE,useConjugacy = TRUE)
Cmodslice <- configureMCMC(mod,print=TRUE,useConjugacy = TRUE, onlySlice = TRUE)
Cmod$removeSamplers(c("reporting","beta","effprop"))
Cmod$addSampler(target = c("reporting","effprop"), type="RW_block")
Cmod$addSampler("beta",type="slice")
Cmod$addMonitors(c("reporting","beta","effprop"))

newMCMC <- buildMCMC(Cmod)
CnewMCMC <- compileNimble(newMCMC,project = mod, resetFunctions = TRUE)

Cnewtime <- system.time(CnewMCMC$run(iterations))
Cnewsample <- as.matrix(CnewMCMC$mvSamples)
effectiveSize(as.mcmc(Cnewsample[,c('beta','effprop','reporting')]))/Cnewtime[3]

```

```{r, "combine everything"}
NimbleCB <- MCMCsuite(code=nimcode,
                   data=nimCBdata,
                   inits=nimCBinits,
                   constants=nimCBcon,
                   MCMCs=c("jags","nimble","nimble_slice","newMCMC"),
                   monitors=c("beta","reporting","effprop"),
                   MCMCdefs = list(
                     newMCMC = quote({
                     Cmod <- configureMCMC(mod,print=FALSE,useConjugacy = TRUE)
                     Cmod$removeSamplers(c("reporting","beta","effprop"))
                     Cmod$addSampler(target = c("reporting","effprop"), type="RW_block")
                     Cmod$addSampler("beta",type="slice")
                     Cmod
                   })),
                   calculateEfficiency=TRUE,
                   niter=iterations,
                   makePlot=FALSE,
                   savePlot=FALSE,
                   setSeed = 5)

print(NimbleCB$timing)
print(NimbleCB$summary)

```

```{r, "conjugates"}
source("nimCB2.R")
mod2 <- nimbleModel(code = nimcode, data=nimCBdata, inits=nimCBinits, constants = nimCBcon,
                    name= "mod2")
Cmod2 <- configureMCMC(mod2,print=TRUE)

NimbleCB2 <- MCMCsuite(code=nimcode,
                   data=nimCBdata,
                   inits=nimCBinits,
                   constants=nimCBcon,
                   MCMCs=c("jags","nimble","nimble_slice"),
                   monitors=c("beta","reporting","effprop"),
                   calculateEfficiency=TRUE,
                   niter=iterations,
                   makePlot=FALSE,
                   savePlot=FALSE,
                   setSeed = 5)
print(NimbleCB$timing)
print(NimbleCB2$timing)
print(NimbleCB2$summary)

mat1 <- as.array(NimbleCB$summary)
mat2 <- as.array(NimbleCB2$summary)
mat1 <- mat1[,"efficiency",]
mat2 <- mat2[,"efficiency",]
a1 <- melt(mat1)
a2 <- melt(mat2)
a1 <- cbind(a1,var3="NoConjugates")
a2 <- cbind(a2,var3="YesConjugates")
dat <- rbind(a1,a2)

ggplot(dat,aes(x=var3, y=value,group=interaction(Var1,Var2))) + geom_line(aes(color=interaction(Var1),linetype=Var2)) + geom_point(aes(color=Var1)) + theme_bw() + ylab('Efficiency') + xlab("speed hack")



```

MCMCsuite can also be used to compare models written in STAN which is described below <br />

<br />
<br />

##### <a name="2.1"> 2.1 "hybrid approach" </a> 

We must rewrite the model so that there are no discrete latent variables. We call this the "hybrid model" <br />
An asside -- Discrete Latent Variables: <br />
An additional asside -- Hamiltonian MCMC: <br />

But before we fit the model in STAN lets explore the hybrid model in NIMBLE <br />

NIMBLE allows us to compare the results of multiple models even if they have different parameterizations 
(e.g. Chain Binomial and the Hybrid Model) <br />

```{r, eval = FALSE}
data$obs <- data$obs + zerohack # Guarnantee that obs remains above 0 (important for the gamma)
data$zerohack <- zerohack

hybridjags <- jags(data = data,
               inits = inits,
               param = params,
               model.file = "hybrid.bug",
               n.iter = 8000,
               n.chains = length(inits))
```

##### <a name="2.1.1"> 2.1.1 Hybrid in JAGS and NIMBLE </a>

```{r, eval = FALSE}
source('nimhybrid.R')
```

```{r, eval = FALSE}
nimhydata <- list(obs = sim$Iobs + zerohack)
nimhycon <- list(numobs = numobs, pop = pop, r0 = r0, zerohack = zerohack)

nimhyinits <- list(I = sim$I + zerohack,
                   effpropS = effpropS,
                   effpropI = effpropI,
                   beta = beta,
                   reporting = reporting,
                   s0 = s0)
```

```{r, eval = FALSE}
nimcb <- MCMCsuite(code = nimcode,
                   data = nimhydata,
                   inits = nimhyinits,
                   constants = nimhycon,
                   MCMCs = c("jags", "nimble"),
                   monitors = c("beta", "reporting", "effpropS", "effpropI"),
                   niter = 10000,
                   makePlot = FALSE,
                   savePlot = FALSE)
```

##### <a name="2.1.2"> 2.1.2 Hybrid in JAGS, NIMBLE and STAN </a>

Run the STAN model <br />

```{r, message=FALSE, verbose = FALSE, warning = FALSE, errors = FALSE, eval = FALSE}
stantime <- system.time (s1 <- stan(file='hybrid.stan', data = data, init = inits,
           pars=c("beta", "reporting", "effpropS", "effpropI", "I"), iter = 8000,
           seed = 1001, chains = length(inits)))
```

Compare all three methods using the hybrid model <br />

```{r, eval = FALSE}
nimhydata <- list(obs = sim$Iobs + zerohack)
nimhycon <- list(numobs = numobs, pop = pop, 
                 r0 = r0, zerohack = zerohack)

nimhyinits <- list(I = sim$I + zerohack,
                   effpropS = effpropS,
                   effpropI = effpropI,
                   beta = beta,
                   reporting = reporting,
                   s0 = s0)
```

```{r, eval = FALSE}
allhybrids <- MCMCsuite(code = nimcode,
                   data = nimhydata,
                   inits = nimhyinits,
                   constants = nimhycon,
                   stan_model = "hybrid.stan",
                   MCMCs = c("jags", "nimble", "stan"),
                   monitors = c("beta", "reporting", "effpropS", "effpropI"),
                   niter = 10000,
                   makePlot = FALSE,
                   savePlot = FALSE)
```

##### <a name="2.2"> 2.2 Finally, compare the Chain Binomial NIMBLE and Hybrid STAN </a>

```{r, eval = FALSE}
nimCBdata <- list(obs = sim$Iobs)
nimCBcon <- list(numobs = numobs, pop = pop, r0 = r0, zerohack = zerohack)

nimCBinits <- list(I = sim$I,
                   effpropS = effpropS,
                   effpropI = effpropI,
                   beta = beta,
                   reporting = reporting,
                   s0 = s0)
```

```{r, eval = FALSE}
nimcb <- MCMCsuite(code = nimcode,
                   data = nimCBdata,
                   inits = nimCBinits,
                   constants = nimCBcon,
                   stan_model = "hybrid.stan",
                   MCMCs = c("jags", "nimble", "stan"),
                   monitors = c("beta", "reporting", "effpropS", "effpropI"),
                   niter = 10000,
                   makePlot = TRUE,
                   savePlot = TRUE)
```

### Part 3

##### <a name="3.1"> 3.1 Expolore more fine-tuned adjustments that can be made in NIMBLE </a>

##### <a name="3.1.1"> 3.1.1 Custom MCMC sampler </a>

The following example is taken directly from pgs 87-89 in the Nimble User Manual <br />

http://r-nimble.org/manuals/ <br />

```{r, eval = FALSE}
## the name of this sampler function, for the purposes of ## adding it to MCMC configurations, will be 'my_RW' my_RW <- nimbleFunction(
## sampler functions must contain 'sampler_BASE'
contains = sampler_BASE,

## sampler functions must have exactly these setup arguments: ## model, mvSaved, target, control
setup = function(model, mvSaved, target, control) {
        ## first, extract the control list elements, which will
        ## dictate the behavior of this sampler.
        ## the setup code will be later processed to determine
        ## all named elements extracted from the control list.
        ## these will become the required elements for any
        ## control list argument to this sampler, unless they also
        ## exist in the NIMBLE system option 'MCMCcontrolDefaultList'.
        ## the random walk proposal standard deviation
        scale <- control$scale
        
        ## determine the list of all dependent nodes,
        ## up to the first layer of stochastic nodes, generally
        ## called 'calcNodes'.  The values, inputs, and logProbs
        ## of these nodes will be retrieved and/or altered
        ## by this algorithm.
calcNodes <- model$getDependencies(target) 
},

## the run function must accept no arguments, execute
## the sampling algorithm, leave the modelValues object
## 'mvSaved' as an exact copy of the updated values in model, ## and have no return value. initially, mvSaved contains
## an exact copy of the values and logProbs in the model.
run = function() {
  
    ## extract the initial model logProb
    model_lp_initial <- getLogProb(model, calcNodes) ## generate a proposal value for target node
    proposal <- rnorm(1, model[[target]], scale)
    ## store this proposed value into the target node.
    ## notice the double assignment operator, `<<-`,
    ## necessary because 'model' is a persistent member
    ## data object of this sampler.
    model[[target]] <<- proposal
    
## calculate target_logProb, propagate the
## proposed value through any deterministic dependents, ## and calculate the logProb for any stochastic
## dependnets. The total (sum) logProb is returned. model_lp_proposed <- calculate(model, calcNodes)
## calculate the log Metropolis-Hastings ratio
    log_MH_ratio <- model_lp_proposed - model_lp_initial
    
## Metropolis-Hastings step: determine whether or ## not to accept the newly proposed value
u <- runif(1, 0, 1)
if (u < exp(log_MH_ratio)) jump <- TRUE
    else jump <- FALSE

## if we accepted the proposal, then store the updated ## values and logProbs from 'model' into 'mvSaved'.
## if the proposal was not accepted, restore the values ## and logProbs from 'mvSaved' back into 'model'. if(jump) copy(from = model, to = mvSaved, row = 1,
nodes = calcNodes, logProb = TRUE) else copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)
},

    ## sampler functions must have a member method 'reset',
    ## which takes no arguments and has no return value.
    ## this function is used to reset the sampler to its
    ## initial state.  since this sampler function maintains
    ## no internal member data variables, reset() needn't
    ## do anything.

    methods = list(reset = function () {} )
)

## now, assume the existence of an R model object 'Rmodel',
## which has a scalar-valued stochastic node 'x'
## create an MCMC configuration with no sampler functions
mcmcspec <- configureMCMC(Rmodel, nodes = NULL)

## add our custom-built random walk sampler on node 'x', ## with a fixed proposal standard deviation = 0.1 mcmcspec$addSampler(target = 'x', type = 'my_RW',
control = list(scale = 0.1))
Rmcmc <- buildMCMC(mcmcspec) ## etc...
        
```

### Part 4

##### <a name="4.1"> NIMBLE extras </a>

##### <a name="4.1.1"> 4.1.1 Mote Carlo Expectation Maximization </a>

Suppose we have a model with missing data (or a layer of latent variables that can be 
treated as missing data) and we would like to maximize the marginal likelihood of the model,
integrating over the missing data. A brute-force method for doing this is MCEM. <br />

Start by building the model <br />

```{r, eval = FALSE}

### Construct the NIMBLE model
CBemout <- nimbleModel(code = nimcode, 
                         name = 'CBemout', 
                         constants = nimCBcon,
                         data = nimCBdata, 
                         inits = nimCBinits)


CBmcem <- buildMCEM(model = CBemout, latentNodes = list("I"), 
                      burnIn = 100, 
                      mcmcControl = list(adaptInterval = 20), 
                      boxConstraints = list( list( c("reporting", "effpropS", "effpropI"), 
                                                   limits = c(0, 1) ) ), 
                        buffer = 1e-6)
# burnIn controls the number of discarded MCMC samples for each "E" step

# The MCEM algorithm allows for box constraints on the nodes that will be optimized,
# specified via the boxConstraints argument. This is highly recommended for nodes that
# have zero density on parts of the real line
```

```{r, eval = FALSE}
CBmcem(maxit = 20, m1 = 150, m2 = 200)
# maxit controls the total number of iterations
# m1 controls the number of MCMC samples for the first half of the iterations
# m2 controls the number of MCMC samples for the second half of the iterations to reduce MCMC error as the algorithm converges
```

Having trouble with "non-finite finite-difference value", so in the meantime here is the example from the user manual <br />
```{r, eval = FALSE}
pumpCode <- nimbleCode({ 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)
```

```{r, eval = FALSE}
box <- list( list(c('alpha','beta'), c(0, Inf)))

pumpMCEM <- buildMCEM(model = pump, latentNodes = 'theta[1:10]',
                       boxConstraints = box)

system.time(pumpMLE <- pumpMCEM())

pumpMLE
```

<br />

##### <a name="4.1.2"> 4.1.2 Particle Filter </a>

Set up the Nimble model <br />

```{r, eval = FALSE}
CBpfout <- nimbleModel(code = nimcode, 
                         name = 'CBpfout', 
                         constants = nimCBcon,
                         data = nimCBdata,
                         check = FALSE)
```

Build the particle filter <br />

```{r, eval = FALSE}
CBpfoutC <- compileNimble(CBpfout)

CBpf <- buildPF(CBpfout, c("I"))
```

```{r, eval = FALSE}
CBpfC <- compileNimble(CBpf, project = CBpfout)
```

Set your parameters <br />

```{r, eval = FALSE}
CBpfoutC$beta = 0.02
CBpfoutC$effpropS = 0.8
CBpfoutC$effpropI = 0.2
CBpfoutC$reporting = 0.5
```

Obtain log-likelihood

```{r, eval = FALSE}
CBpfC$run(m = 5000)
```

Currently relatively useless as is... <br />

Use this framework to construct your own updater <br />

<br />

##### Part 5:  

<a name="5.1"> NIMBLE Notes </a>

**Truncation of distributions** <br />
   • x ∼ N(0, sd = 10) T(0, a)$, or <br />
   • x ∼ T(dnorm(0, sd = 10), 0, a), <br />

   • mu1 ~ dnorm(0, 1) <br />
   • mu2 ~ dnorm(0, 1) <br />
   • constraint_data ~ dconstraint( mu1 + mu2 > 0 ) <br />
   
**Lifted Nodes** <br />
   • The use of link functions causes new nodes to be introduced <br />
   • When distribution parameters are expressions, NIMBLE creates a new deterministic node that contains the expression for a given parameter <br />
   
**logProb** <br />
    • For each variable that contains at least one stochastic node, NIMBLE generates a model variable with the prefix “logProb” <br />
    • Can be retrieved with getLogProb
    
**Choice of Samplers** <br /> 
    1. If the node has no stochastic dependents, a predictive end sampler is assigned. The end sampling algorithm merely calls simulate on the particular node. <br /> 
    2. The node is checked for presence of a conjugate relationship between its prior distribution and the distributions of its stochastic dependents. If it is determined to be in a conjugate relationship, then the corresponding conjugate (Gibbs) sampler is assigned. <br /> 
    3. If the node is discrete-valued, then a slice sampler is assigned [5]. <br /> 
    4. If the node follows a multivariate distribution, then a RW block sampler is assigned for all elements. This is a Metropolis-Hastings adaptive random-walk sampler with a multivariate normal proposal [6]. <br /> 
    5. If none of the above criteria are satisfied, then a RW sampler is assigned. This is a Metropolis-Hastings adaptive random-walk sampler with a univariate normal proposal distribution. <br /> 
    
**Missing Values** <br /> 

```{r, eval = FALSE}
CBMiss <- CBout$newModel()
CBMiss$resetData()
CBDataNew <- nimCBdata
CBDataNew$obs[c(1, 3)] <- NA
CBMiss$setData(CBDataNew)
CBMissSpec <- configureMCMC(CBMiss)
CBMissSpec$addMonitors(c("obs", "beta", "effpropS", "effpropI", "reporting"))
CBMissMCMC <- buildMCMC(CBMissSpec)
CBCobj <- compileNimble(CBMiss, CBMissMCMC)
niter <- 1000
set.seed(0)
CBCobj$CBMissMCMC$run(niter)
samples <- as.matrix(CBCobj$CBMissMCMC$mvSamples)
out <- rbind(head(samples), tail(samples))
out
```

**Multiple instances of the same model** <br />
Sometimes it is useful to have more than one copy of the same model. For example, nimbleFunctions are often bound to a particular model as a result of setup code.
```{r, eval = FALSE}
simpleCode <- nimbleCode({
for(i in 1:N) x[i] ~ dnorm(0, 1) })
## Return the model definition only, not a built model
simpleModelDefinition <- nimbleModel(simpleCode, constants = list(N = 10),
                                     returnDef = TRUE, check = FALSE)

## Make one instance of the model
simpleModelCopy1 <- simpleModelDefinition$newModel(check = FALSE) ## Make another instance from the same definition
simpleModelCopy2 <- simpleModelDefinition$newModel(check = FALSE) ## Ask simpleModelCopy2 for another copy of itself
simpleModelCopy3 <- simpleModelCopy2$newModel(check = FALSE)
```