# Gp Nimble Setup

library("nimble")
library("ggplot2")
library("tidyr")
library("dplyr")

## dmnorm example

from Daniel, who notes we need to be explicit about vector/matrix sizes on both LHS and RHS,

code <- nimbleCode({
mu ~ dunif(-100, 100)
sigma ~ dunif(0, 100)
rho ~ dunif(20, 100)
muVec[1:N] <- mu * onesVector[1:N]
Cov[1:N, 1:N] <- sigma^2 * exp(-dist[1:N, 1:N]/rho)
g[1:N] ~ dmnorm(muVec[1:N], cov = Cov[1:N, 1:N])
for (i in 1:N) {
y[i] ~ dpois(exp(g[i]))
}
})
constants <- list(N = 148,
onesVector = rep(1, 148),
dist = matrix(1, nrow=148, ncol=148))
y = rpois(148, 1) # vector of 148 observations (counts)
data <- list(y=y)

inits <- list(mu = 0,
sigma = 5,
rho = 60,
g = rep(0, 148))
m <- nimbleModel(code = code, constants = constants, data = data, inits = inits)

## GP Model, original specification

This model definition doesn’t work, since we do not define sizes explicitly

code <- nimbleCode({
l ~ dunif(0,1e4)
sigma.n ~ dunif(0, 1e4)
Sigma[1:N, 1:N] <- exp(-0.5 * diff / l ^ 2) + sigma.n ^ 2 * I
y[1:N] ~ dmnorm(Mu, cov = Sigma)
})
obs <- data.frame(x = c(-4, -3, -1,  0,  2),
y = c(-2,  0,  1,  2, -1))

N <- length(obs$x) diff <- outer(obs$x, obs$x, function(xi, xj) (xi-xj)^2) constants = list(N = N, x = obs$x, Mu = rep(0,N), diff = diff, I = diag(1,N))
inits <- list(l = 1, sigma.n = 10)
data <- data.frame(y = obs$y) m <- nimbleModel(code = code, constants = constants, data = data, inits = inits) Error in x$parentNodeExprs[[iV]][[iDim + 2]]: object of type 'symbol' is not subsettable

## GP Model, v2

This model definition still doesn’t work, even though it appears to be more explicit

code <- nimbleCode({

l ~ dunif(0,100)
sigma.n ~ dunif(0,100)
SE <- function(xi,xj, l) exp(-0.5 * (xi - xj) ^ 2 / l ^ 2)
Sigma[1:N, 1:N] <- outer(x[1:N], x[1:N], SE, l) + sigma.n ^ 2 * I[1:N, 1:N]
y[1:N] ~ dmnorm(Mu[1:N], cov = Sigma[1:N, 1:N])

})
obs <- data.frame(x = c(-4, -3, -1,  0,  2),
y = c(-2,  0,  1,  2, -1))
N <- length(obs$x) constants = list(N = N, x = obs$x, Mu = rep(0,N), I = diag(1,N))
inits <- list(l = 1, sigma.n = 10)
data <- data.frame(y = obs$y) nimbleModel(code = code, constants = constants, data = data, inits = inits) Error in replaceConstantsRecurse(x, constEnv, constNames): Error, hit end ## GP Model A successful specification looks like this: code <- nimbleCode({ l ~ dunif(0,1e5) # dgamma(5, 5) sigma.n ~ dunif(0, 1e5) # dgamma(5, 5) Sigma[1:N, 1:N] <- exp(-0.5 * diff[1:N, 1:N] / l ^ 2) + sigma.n ^ 2 * I[1:N, 1:N] y[1:N] ~ dmnorm(Mu[1:N], cov = Sigma[1:N, 1:N]) }) obs <- data.frame(x = c(-4, -3, -1, 0, 2), y = c(-2, 0, 1, 2, -1)) N <- length(obs$x)
diff <- outer(obs$x, obs$x, function(xi, xj) (xi-xj)^2)
constants = list(N = N, x = obs$x, Mu = rep(0,N), diff = diff, I = diag(1,N)) inits <- list(l = 1, sigma.n = 10) data <- data.frame(y = obs$y)
m <- nimbleModel(code = code, constants = constants, data = data, inits = inits)

Define our mcmc procedure in Nimble

my_mcmc <- function(code, constants, data, inits, n_iter=1e5, thin = 1e2){
Rmodel <- nimbleModel(code = code, constants = constants, data = data, inits = inits)
Cmodel <- compileNimble(Rmodel)
mcmcspec <- configureMCMC(Rmodel, print=FALSE,thin=thin)
Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
Cmcmc$run(n_iter) samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
#samples <- samples[,1:(length(inits)-1)]
gather(samples)
}
df <- my_mcmc(code = code, constants = constants, data = data, inits = inits, n = 1e6)

#### Summary statistics

summarise(group_by(df, key), mean=mean(value), std=sqrt(var(value)))
Source: local data frame [2 x 3]

key         mean          std
1       l 49845.347322 29128.275365
2 sigma.n     2.188315     1.372436

#### Traces

ggplot(sample_n(df,500)) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')

#### Posteriors

ggplot(df) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')

## Time-series example

f <- function (x, p){ x * exp(p[1] * (1 - x / p[2]) * (x - p[3]) / p[2]) }
p <- c(2, 8, 5)
z_g <- function() rlnorm(1, 0, 0.05)
T <- 40

set.seed(1234)
x <- numeric(T)
x[1] <- 5.5

for(t in 1:(T-1))
x[t+1] = z_g() * f(x[t], p=p)

qplot(seq_along(x), x)
data <- data.frame(x=x)

#### model setup

obs <- data.frame(x=x[1:(T-1)], y = x[2:T])
N <- length(obs$x) diff <- outer(obs$x, obs$x, function(xi, xj) (xi-xj)^2) constants = list(N = N, x = obs$x, Mu = rep(0,N), diff = diff, I = diag(1,N))
inits <- list(l = 1, sigma.n = 10)
data <- data.frame(y = obs\$y)

#### Estimate GP bm MCMC with Nimble

df <- my_mcmc(code = code, constants = constants, data = data, inits = inits, n = 1e6)

#### Summary statistics

summarise(group_by(df, key), mean=mean(value), std=sqrt(var(value)))
Source: local data frame [2 x 3]

key      mean        std
1       l 3.8059775 0.60843921
2 sigma.n 0.3641306 0.04565619

#### Traces

ggplot(sample_n(df,500)) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')

#### Posteriors

ggplot(df) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')