Lsn Nimble

some sample data

library(sde)
library(nimble)
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=20, N=100))

sigma.x not provided, attempting symbolic derivation.

plot(data)

LSN version

Test case: Set prior for m \(\approx 0\):

lsn <- modelCode({
   theta ~ dunif(1e-10, 100.0)
   sigma_x ~ dunif(1e-10, 100.0)
   sigma_y ~ dunif(1e-10, 100.0)
       m ~ dunif(-1e2, 1e2)
    x[1] ~ dunif(0, 100)
    y[1] ~ dunif(0, 100) 

  for(i in 1:(N-1)){
    mu_x[i] <- x[i] + y[i] * (theta - x[i]) 
    x[i+1] ~ dnorm(mu_x[i], sd = sigma_x) 
    mu_y[i] <- y[i] + m * t[i]
    y[i+1] ~ dnorm(mu_y[i], sd = sigma_y) 
  }
})

Constants in the model definition are the length of the dataset, \(N\) and the time points of the sample. Note we’ve made time explicit, we’ll assume uniform spacing here.

constants <- list(N = length(data$x), t = 1:length(data$x))

Initial values for the parameters

inits <- list(theta = 6, m = 0, sigma_x = 1, sigma_y = 1, y = rep(1,constants$N))

and here we go as before:

Rmodel <- nimbleModel(code = lsn, 
                      constants = constants, 
                      data = data, 
                      inits = inits)
Cmodel <- compileNimble(Rmodel)

mcmcspec <- MCMCspec(Rmodel, print=TRUE,thin=1e2)

[1] RW sampler;   targetNode: theta,  adaptive: TRUE,  adaptInterval: 200,  scale: 1
[2] RW sampler;   targetNode: sigma_x,  adaptive: TRUE,  adaptInterval: 200,  scale: 1
[3] RW sampler;   targetNode: sigma_y,  adaptive: TRUE,  adaptInterval: 200,  scale: 1
[4] RW sampler;   targetNode: m,  adaptive: TRUE,  adaptInterval: 200,  scale: 1
[5] RW sampler;   targetNode: y[1],  adaptive: TRUE,  adaptInterval: 200,  scale: 1
[6] conjugate_dnorm sampler;   targetNode: y[2],  dependents_dnorm: x[3], y[3]
...

Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)

Cmcmc(1e4)

NULL

and examine results

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

[1] 100 206

samples <- samples[,1:4]

mean(samples$theta)

[1] 10.11174

mean(samples$m)

[1] -1.88765e-05

mean(samples$sigma_x)

[1] 0.385018

plot(samples[ , 'm'], type = 'l', xlab = 'iteration', ylab = 'm')
plot(samples[ , 'sigma_x'], type = 'l', xlab = 'iteration', ylab = expression(sigma[x]))
plot(samples[ , 'sigma_y'], type = 'l', xlab = 'iteration', ylab = expression(sigma[y]))
plot(samples[ , 'theta'], type = 'l', xlab = 'iteration', ylab = expression(theta))

hist(samples[, 'm'], xlab = 'm')
hist(samples[, 'sigma_x'], xlab = expression(sigma[x]))
hist(samples[, 'sigma_y'], xlab = expression(sigma[y]))
hist(samples[, 'theta'], xlab = expression(theta))