Thursday_NIMBLE.Rmd

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

```{r setup}
library(knitr)
opts_knit$set(root.dir = "/Users/Morgan/Documents/School Work/McMaster/NIMBLE_Presentation/")
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")
```

<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 = 3)
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

```{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(1, 2), mai = c(.6, .5, .1, .2))
  plot(samples[ , 'beta'], type = 'l', xlab = 'iteration')
  plot(samples[ , 'effpropS'], type = 'l', xlab = 'iteration')
  plot(samples[ , 'effpropI'], type = 'l', xlab = 'iteration')
  plot(samples[ , 'reporting'], type = 'l', xlab = 'iteration')
  par(mfrow = c(1, 1))
  plot(samples[ , 'effpropS'], samples[ , 'effpropI'], pch = 20)
```

Save these points for later

```{r}
  def_effS <- samples[ , 'effpropS']
  def_effI <- samples[ , 'effpropI']
```

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(30)
```

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

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

```{r}
CBoutMCMC$run(30000)
samplesNew <- as.matrix(CBoutMCMC$mvSamples)
```

Check for an imporvement

```{r}
  par(mfrow = c(2,2))
  acf(samplesNew[, "effpropS"])
  acf(samplesNew[, "effpropI"])
  plot(samplesNew[ , 'effpropS'], type = 'l', xlab = 'iteration')
  plot(samplesNew[ , 'effpropI'], type = 'l', xlab = 'iteration')
  par(mfrow = c(1,1))
  plot(samplesNew[ , 'effpropS'], samplesNew[ , 'effpropI'], pch = 20)
  points(def_effS, def_effI, pch = 20, col = "blue")
```

Well that didn't do anything... <br />

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

```{r}
CBout <- nimbleModel(code = nimcode, 
                         name = 'CBout', 
                         constants = nimCBcon,
                         data = nimCBdata, 
                         inits = nimCBinits)
  
CBoutC <- compileNimble(CBout)
  
CBoutSpec <- configureMCMC(CBout, print = TRUE)

CBoutSpec$addMonitors(c("beta", "effpropS", "effpropI", "reporting"))

CBoutSpec$setThin(20)
```

```{r}
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)
```

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

```{r}
nimtimevec[7] <- system.time(CBoutMCMC$run(30000))[3]
samplesblock <- as.matrix(CBoutMCMC$mvSamples)
nim_eff <- effectiveSize(as.mcmc.list(as.mcmc(samplesblock))) / nimtimevec[7]
nim_eff
```

```{r}
  par(mfrow = c(1,1))
  plot(samplesblock[ , 'effpropS'], samplesblock[ , 'effpropI'], pch = 20)
```

<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}
nimcb <- MCMCsuite(code = nimcode,
                   data = nimCBdata,
                   inits = nimCBinits,
                   constants = nimCBcon,
                   MCMCs = c("jags", "nimble"),
                   monitors = c("beta", "reporting", "effpropS", "effpropI"),
                   niter = 12000,
                   calculateEfficiency = TRUE,
                   makePlot = FALSE,
                   savePlot = FALSE)
  
  nimcb$summary
  
  nimcb$efficiency
  
  nimcb$timing
```

MCMCsuite allows for convenient specification of MCMC runs <br />

For example, we may want to run a bit longer, introduce some burn-in, thinning and also compare a adaptive block sampler <br />

```{r}
nimcb <- MCMCsuite(code = nimcode,
                   data = nimCBdata,
                   inits = nimCBinits,
                   constants = nimCBcon,
                   MCMCs = c("jags", "nimble", "autoBlock", "nimble_RW"),
                   monitors = c("beta", "reporting", "effpropS", "effpropI"),
                   niter = 12000,
                   burnin = 1000,
                #   thin = 10,
                   calculateEfficiency = TRUE,
                   makePlot = FALSE,
                   savePlot = FALSE)
  
  nimcb$summary
  
  nimcb$efficiency
  
  nimcb$timing
```

MCMCsuite can also be used to compare models written in STAN or... which we will see next week <br />

<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}
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}
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)
```